fix: handle small negative rates in splitter due to floating precision issues (#1305)

tests: add test of splitter

chore: add docstrings

Refs.
equinor/ecalc-internal#1372

Author: olelod

docs/drafts/next.draft.md

@@ -12,6 +12,8 @@ sidebar_position: -58
 
 ## Bug Fixes
 
+- Fixed an issue where floating point precision errors could lead to negative mass rates within multiple streams and pressures compressor trains under certain conditions.
+
 ## Breaking changes
 
 ### CLI

src/libecalc/domain/process/compressor/core/train/compressor_train_common_shaft_multiple_streams_and_pressures.py

@@ -480,15 +480,15 @@ def calculate_compressor_train(
         for stage_number, stage in enumerate(self.stages):
             stage_inlet_stream = previous_stage_outlet_stream
 
-            additional_rates_to_splitter = [
+            rates_out_of_splitter = [
                 constraints.stream_rates[stream_number]
                 for stream_number in self.outlet_stream_connected_to_stage.get(stage_number, [])
             ]
-            additional_streams_to_mixer = []
+            streams_in_to_mixer = []
             for stream_number in self.inlet_stream_connected_to_stage.get(stage_number, []):
                 if stream_number > 0:
                     if inlet_stream_counter < len(fluid_streams):
-                        additional_streams_to_mixer.append(
+                        streams_in_to_mixer.append(
                             fluid_streams[inlet_stream_counter].create_stream_with_new_conditions(
                                 conditions=ProcessConditions(
                                     pressure_bara=stage_inlet_stream.pressure_bara,
@@ -501,8 +501,8 @@ def calculate_compressor_train(
             stage_results.append(
                 stage.evaluate(
                     inlet_stream_stage=stage_inlet_stream,
-                    additional_rates_to_splitter=additional_rates_to_splitter,
-                    additional_streams_to_mixer=additional_streams_to_mixer,
+                    rates_out_of_splitter=rates_out_of_splitter,
+                    streams_in_to_mixer=streams_in_to_mixer,
                     speed=self.shaft.get_speed(),
                     asv_rate_fraction=asv_rate_fraction,
                     asv_additional_mass_rate=asv_additional_mass_rate,

src/libecalc/domain/process/compressor/core/train/stage.py

@@ -141,8 +141,8 @@ def evaluate(
         self,
         inlet_stream_stage: FluidStream,
         speed: float,
-        additional_rates_to_splitter: list[float] | None = None,
-        additional_streams_to_mixer: list[FluidStream] | None = None,
+        rates_out_of_splitter: list[float] | None = None,
+        streams_in_to_mixer: list[FluidStream] | None = None,
         asv_rate_fraction: float | None = 0.0,
         asv_additional_mass_rate: float | None = 0.0,
     ) -> CompressorTrainStageResultSingleTimeStep:
@@ -153,8 +153,8 @@ def evaluate(
             inlet_stream_stage (FluidStream): The conditions of the inlet fluid stream. If there are several inlet streams,
                 the first one is the stage inlet stream, the others enter the stage at the Mixer.
             speed (float): The speed of the shaft driving the compressor
-            additional_rates_to_splitter (list[float] | None, optional): Additional rates to the Splitter if defined.
-            additional_streams_to_mixer (list[FluidStream] | None, optional): Additional streams to the Mixer if defined.
+            rates_out_of_splitter (list[float] | None, optional): Additional rates to the Splitter if defined.
+            streams_in_to_mixer (list[FluidStream] | None, optional): Additional streams to the Mixer if defined.
             asv_rate_fraction (float | None, optional): Fraction of the available capacity of the compressor to fill
                 using some kind of pressure control (on the interval [0,1]). Defaults to 0.0.
             asv_additional_mass_rate (float | None, optional): Additional recirculated mass rate due to
@@ -165,22 +165,20 @@ def evaluate(
         """
         # First the stream passes through the Splitter (if defined)
         if self.splitter is not None:
-            if additional_rates_to_splitter is None:
-                raise IllegalStateException("additional_rates_to_splitter cannot be None when a splitter is defined")
+            self.splitter.rates_out_of_splitter = rates_out_of_splitter
             inlet_stream_after_splitter = self.split(
                 inlet_stream_stage=inlet_stream_stage,
-                additional_rates_to_splitter=additional_rates_to_splitter,
             )
         else:
             inlet_stream_after_splitter = inlet_stream_stage
 
         # Then the stream passes through the Mixer     (if defined)
         if self.mixer is not None:
-            if additional_streams_to_mixer is None:
-                raise IllegalStateException("additional_streams_to_mixer cannot be None when a mixer is defined")
+            if streams_in_to_mixer is None:
+                raise IllegalStateException("streams_in_to_mixer cannot be None when a mixer is defined")
             inlet_stream_after_mixer = self.mix(
                 inlet_stream_stage=inlet_stream_after_splitter,
-                additional_streams_to_mixer=additional_streams_to_mixer,
+                streams_in_to_mixer=streams_in_to_mixer,
             )
         else:
             inlet_stream_after_mixer = inlet_stream_after_splitter
@@ -257,38 +255,26 @@ def evaluate(
     def split(
         self,
         inlet_stream_stage: FluidStream,
-        additional_rates_to_splitter: list[float],
     ) -> FluidStream:
         """Split the inlet stream into many streams. One stream goes to the compressor stage. The other(s) are taken out.
         In the future, the additional streams could be used for other purposes, but today they are just dropped completely.
 
         Args:
             inlet_stream_stage (FluidStream): The inlet stream for the stage.
-            additional_rates_to_splitter (list[float]): Additional rates to split from the inlet stream.
 
         Returns:
             FluidStream: The stream going to the compressor stage.
         """
         assert self.splitter is not None
-        assert additional_rates_to_splitter is not None
-        if self.splitter.number_of_outputs != len(additional_rates_to_splitter) + 1:
-            raise IllegalStateException(
-                f"Number of additional rates to Splitter ({len(additional_rates_to_splitter)}) "
-                f"does not match number of Splitter outputs ({self.splitter.number_of_outputs})."
-            )
-        all_rates_to_splitter = additional_rates_to_splitter + [
-            inlet_stream_stage.standard_rate - sum(additional_rates_to_splitter)
-        ]
         split_streams = self.splitter.split_stream(
             stream=inlet_stream_stage,
-            split_fractions=all_rates_to_splitter,
         )
         return split_streams[-1]  # The last stream goes to the compressor stage
 
     def mix(
         self,
         inlet_stream_stage: FluidStream,
-        additional_streams_to_mixer: list[FluidStream],
+        streams_in_to_mixer: list[FluidStream],
         prefer_first_stream: bool = True,
     ) -> FluidStream:
         """Mix the inlet stream with additional streams.
@@ -300,22 +286,22 @@ def mix(
 
         Args:
             inlet_stream_stage (FluidStream): The inlet stream for the stage.
-            additional_streams_to_mixer (list[FluidStream]): Additional streams to mix with the inlet stream.
+            streams_in_to_mixer (list[FluidStream]): Additional streams to mix with the inlet stream.
             prefer_first_stream (bool): Whether to prefer the properties of the first stream when mixing streams
                                         with zero mass flow. Defaults to True. (Which fluid to recirculate)
 
         Returns:
             FluidStream: The mixed stream.
         """
         assert self.mixer is not None
-        assert additional_streams_to_mixer is not None
-        if self.mixer.number_of_inputs != len(additional_streams_to_mixer) + 1:
+        assert streams_in_to_mixer is not None
+        if self.mixer.number_of_inputs != len(streams_in_to_mixer) + 1:
             raise IllegalStateException(
-                f"Number of additional rates to Splitter ({len(additional_streams_to_mixer)}) "
+                f"Number of additional rates to Splitter ({len(streams_in_to_mixer)}) "
                 f"does not match number of Splitter outputs ({self.splitter.number_of_inputs})."
             )
 
-        all_streams_to_mixer = [inlet_stream_stage] + additional_streams_to_mixer
+        all_streams_to_mixer = [inlet_stream_stage] + streams_in_to_mixer
         if sum(s.mass_rate_kg_per_h for s in all_streams_to_mixer) == 0:
             if prefer_first_stream:
                 return inlet_stream_stage

src/libecalc/domain/process/compressor/core/train/utils/common.py

@@ -4,6 +4,7 @@
 from libecalc.domain.process.compressor.core.train.utils.numeric_methods import DampState, adaptive_pressure_update
 from libecalc.domain.process.value_objects.fluid_stream import FluidStream
 
+FLOATING_POINT_PRECISION = 1e-6
 EPSILON = 1e-5
 PRESSURE_CALCULATION_TOLERANCE = 1e-3
 POWER_CALCULATION_TOLERANCE = 1e-3

src/libecalc/domain/process/entities/process_units/splitter/splitter.py

@@ -1,33 +1,92 @@
+from libecalc.common.logger import logger
 from libecalc.domain.process.value_objects.fluid_stream import FluidStream
 
 
 class Splitter:
-    def __init__(self, number_of_outputs: int):
-        self.number_of_outputs = number_of_outputs
+    """
+    A class representing a splitter that divides a fluid stream into multiple output streams.
+
+    Attributes:
+        number_of_outputs (int): The number of output streams.
+        _rates_out_of_splitter (list[float]): The flow rates for the first (number_of_outputs - 1) outputs.
+    """
 
-    def split_stream(self, stream: FluidStream, split_fractions: list[float]) -> list[FluidStream]:
+    def __init__(self, number_of_outputs: int, rates_out_of_splitter: list[float] | None = None):
         """
-        Splits the fluid stream into streams based on the split fractions.
+        Initializes the Splitter instance.
 
         Args:
-            stream (FluidStream): The fluid stream to be split.
-            split_fractions (Sequence[float]): The fractions of the stream to go to the different output streams (0.0 to 1.0).
+            number_of_outputs (int): The number of output streams. Must be larger than 0.
+            rates_out_of_splitter (list[float] | None): The flow rates for the first (number_of_outputs - 1) outputs.
+                If None, initializes with zero flow rates.
+
+        Raises:
+            ValueError: If the length of rates_out_of_splitter does not match (number_of_outputs - 1).
+        """
+        self.number_of_outputs = number_of_outputs
+        assert isinstance(self.number_of_outputs, int) and self.number_of_outputs > 0
+
+        self._rates_out_of_splitter = (
+            [0.0] * (number_of_outputs - 1) if rates_out_of_splitter is None else rates_out_of_splitter
+        )
+        assert len(self._rates_out_of_splitter) == self.number_of_outputs - 1
+
+    @property
+    def rates_out_of_splitter(self) -> list[float]:
+        """
+        Gets the flow rates for the first (number_of_outputs - 1) outputs.
 
         Returns:
-            list[FluidStream]: A list containing the resulting FluidStreams.
+            list[float]: The flow rates for the outputs.
+        """
+        return self._rates_out_of_splitter
+
+    @rates_out_of_splitter.setter
+    def rates_out_of_splitter(self, rates: list[float] | None):
+        """
+        Sets the flow rates for the first (number_of_outputs - 1) outputs.
+
+        Args:
+            rates (list[float] | None): The new flow rates. If None, initializes with zero flow rates.
+
+        Raises:
+            ValueError: If the length of rates does not match (number_of_outputs - 1).
         """
-        # make sure number of split fractions matches number of outputs
-        if len(split_fractions) != self.number_of_outputs:
-            raise ValueError("Number of split fractions must match number of outputs.")
+        self._rates_out_of_splitter = [0.0] * (self.number_of_outputs - 1) if rates is None else rates
+        assert len(self._rates_out_of_splitter) == self.number_of_outputs - 1
+
+    def split_stream(self, stream: FluidStream) -> list[FluidStream]:
+        """
+        Splits the input fluid stream into multiple output streams.
+
+        The flow rates for the first (number_of_outputs - 1) outputs are taken from rates_out_of_splitter.
+        The flow rate for the last output is calculated to ensure mass balance. If the sum of rates_out_of_splitter
+        exceeds the input stream rate (checked against FLOATING_POINT_PRECISION), an exception is raised.
 
-        # normalize split fractions
-        total = sum(split_fractions)
-        normalized_fractions = [f / total for f in split_fractions]
+        Args:
+            stream (FluidStream): The fluid stream to be split.
+
+        Returns:
+            list[FluidStream]: A list of FluidStream objects representing the output streams.
 
+        Raises:
+            IllegalStateException: If the sum of rates_out_of_splitter exceeds the input stream rate.
+        """
+        last_rate = stream.standard_rate - sum(self.rates_out_of_splitter)
+        if last_rate < 0:
+            logger.warning(
+                f"Sum of output rates from splitter slightly exceeds input stream rate, probably due to floating point precision. "
+                f"Input stream rate: {stream.standard_rate}, "
+                f"Sum of output rates: {sum(self.rates_out_of_splitter)}. "
+                f"Correcting last output rate to zero."
+            )
+            last_rate = 0.0
+        all_rates_out_of_splitter = self.rates_out_of_splitter + [last_rate]
+        split_fractions = [rate / sum(all_rates_out_of_splitter) for rate in all_rates_out_of_splitter]
         return [
             FluidStream(
                 thermo_system=stream.thermo_system,
                 mass_rate_kg_per_h=stream.mass_rate_kg_per_h * split_fraction,
             )
-            for split_fraction in normalized_fractions
+            for split_fraction in split_fractions
         ]

tests/libecalc/domain/process/entities/process_units/test_splitter.py

@@ -0,0 +1,40 @@
+from ecalc_neqsim_wrapper.thermo import STANDARD_PRESSURE_BARA, STANDARD_TEMPERATURE_KELVIN
+from libecalc.domain.process.entities.process_units.splitter.splitter import Splitter
+from libecalc.domain.process.value_objects.fluid_stream import FluidComposition
+from libecalc.domain.process.value_objects.fluid_stream.fluid_model import EoSModel, FluidModel
+from libecalc.infrastructure.neqsim_fluid_provider.neqsim_fluid_factory import NeqSimFluidFactory
+
+
+def test_splitter():
+    composition = FluidComposition(
+        nitrogen=3,
+        CO2=1,
+        methane=62,
+        ethane=15,
+        propane=13,
+        i_butane=1,
+        n_butane=2,
+        i_pentane=1,
+        n_pentane=1,
+        n_hexane=1,
+        water=25,
+    )
+    inlet_stream = NeqSimFluidFactory(
+        FluidModel(
+            eos_model=EoSModel.SRK,
+            composition=composition,
+        )
+    ).create_stream_from_standard_rate(
+        pressure_bara=STANDARD_PRESSURE_BARA,
+        temperature_kelvin=STANDARD_TEMPERATURE_KELVIN,
+        standard_rate_m3_per_day=100000,
+    )
+    splitter = Splitter(
+        number_of_outputs=2,
+        rates_out_of_splitter=[50000],
+    )
+    outlet_streams = splitter.split_stream(inlet_stream)
+
+    assert inlet_stream.standard_rate == 100000
+    assert outlet_streams[0].standard_rate == 50000
+    assert outlet_streams[1].standard_rate == 50000