diff --git a/docs/drafts/next.draft.md b/docs/drafts/next.draft.md index 88ef8ad32b..ba37489cba 100644 --- a/docs/drafts/next.draft.md +++ b/docs/drafts/next.draft.md @@ -12,6 +12,7 @@ sidebar_position: -61 STP: "flare" column has been added to STP Export - for `FIXED` installations only. ## Bug Fixes +Fixed: `SpeedSolver` now correctly handles EOS/PHflash failures near the dense/supercritical boundary instead of propagating unhandled exceptions. - Hardened compressor PH flash handling so invalid thermodynamic states are no longer used in compressor outlet calculations. diff --git a/src/libecalc/process/process_pipeline/process_error.py b/src/libecalc/process/process_pipeline/process_error.py index 9ac05b31ab..7c294a8157 100644 --- a/src/libecalc/process/process_pipeline/process_error.py +++ b/src/libecalc/process/process_pipeline/process_error.py @@ -1,7 +1,25 @@ +from dataclasses import dataclass + from libecalc.common.errors.exceptions import EcalcError from libecalc.process.process_pipeline.process_unit import ProcessUnitId +@dataclass(frozen=True) +class CompressorOperatingPoint: + """Thermodynamic and mechanical state describing what a compressor was asked to do. + + Together these fields are sufficient to reproduce the operating point in isolation + (build the inlet stream, set the speed, request the head/efficiency). + """ + + inlet_pressure_bara: float + inlet_temperature_kelvin: float + actual_rate_m3_per_hour: float + polytropic_head_joule_per_kg: float + polytropic_efficiency: float + speed: float + + class ProcessError(EcalcError): def __init__(self, reason: str | None = None): self.reason = reason @@ -27,6 +45,29 @@ def __init__( super().__init__(reason) +class OutletFluidNotAchievableError(ProcessError): + """Raised when the compressor cannot produce a thermodynamically valid outlet stream. + + This typically means the EOS / PH flash rejected the requested outlet state + (e.g. NaN/inf properties, enthalpy did not converge, or NeqSim raised an + exception). Unlike RateTooHighError the compressor was operating inside its + chart capacity — the fluid itself is the problem. + + All fields describe the exact operating point at failure, enough to reproduce + the issue in isolation. + """ + + def __init__( + self, + process_unit_id: ProcessUnitId, + unachievable_operating_point: CompressorOperatingPoint, + reason: str = "Outlet fluid state is not achievable.", + ): + self.process_unit_id = process_unit_id + self.unachievable_operating_point = unachievable_operating_point + super().__init__(reason) + + class RateTooHighError(OutsideCapacityError): def __init__( self, diff --git a/src/libecalc/process/process_solver/search_strategies.py b/src/libecalc/process/process_solver/search_strategies.py index 0b198623e0..4e33de6a35 100644 --- a/src/libecalc/process/process_solver/search_strategies.py +++ b/src/libecalc/process/process_solver/search_strategies.py @@ -1,5 +1,6 @@ import abc from collections.abc import Callable +from typing import NamedTuple from scipy.optimize import root_scalar @@ -9,6 +10,17 @@ CONVERGENCE_TOLERANCE = 1e-5 +class BinarySearchResult(NamedTuple): + higher: bool + accepted: bool + + +ACCEPT_AND_GO_HIGHER = BinarySearchResult(higher=True, accepted=True) +REJECT_AND_GO_LOWER = BinarySearchResult(higher=False, accepted=False) +ACCEPT_AND_GO_LOWER = BinarySearchResult(higher=False, accepted=True) +REJECT_AND_GO_HIGHER = BinarySearchResult(higher=True, accepted=False) + + class DidNotConvergeError(EcalcError): def __init__( self, diff --git a/src/libecalc/process/process_solver/solver.py b/src/libecalc/process/process_solver/solver.py index 5ac9caa614..b61c7afef5 100644 --- a/src/libecalc/process/process_solver/solver.py +++ b/src/libecalc/process/process_solver/solver.py @@ -8,6 +8,7 @@ from typing import Self, TypeVar from libecalc.process.fluid_stream.fluid_stream import FluidStream +from libecalc.process.process_pipeline.process_error import CompressorOperatingPoint from libecalc.process.process_pipeline.process_pipeline import ProcessPipelineId from libecalc.process.process_pipeline.process_unit import ProcessUnitId from libecalc.process.process_solver.configuration import ( @@ -25,6 +26,7 @@ class SolverFailureStatus(StrEnum): BELOW_MINIMUM_FLOW_RATE = "BELOW_MINIMUM_FLOW_RATE" MAXIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_BELOW_TARGET = "MAXIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_BELOW_TARGET" MINIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_ABOVE_TARGET = "MINIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_ABOVE_TARGET" + OUTLET_FLUID_NOT_ACHIEVABLE = "OUTLET_FLUID_NOT_ACHIEVABLE" @dataclass @@ -46,7 +48,14 @@ def with_source_id(self, source_id: ProcessPipelineId) -> Self: return dataclasses.replace(self, source_id=source_id) -SolverFailureEvent = OutsideCapacityEvent | TargetNotAchievableEvent +@dataclass +class OutletFluidNotAchievableEvent: + status: SolverFailureStatus + source_id: ProcessUnitId + unachievable_operating_point: CompressorOperatingPoint + + +SolverFailureEvent = OutsideCapacityEvent | TargetNotAchievableEvent | OutletFluidNotAchievableEvent @dataclass(frozen=True) diff --git a/src/libecalc/process/process_solver/solvers/speed_solver.py b/src/libecalc/process/process_solver/solvers/speed_solver.py index b48e7e12c4..f3e23b7246 100644 --- a/src/libecalc/process/process_solver/solvers/speed_solver.py +++ b/src/libecalc/process/process_solver/solvers/speed_solver.py @@ -2,11 +2,24 @@ from collections.abc import Callable from libecalc.process.fluid_stream.fluid_stream import FluidStream -from libecalc.process.process_pipeline.process_error import RateTooHighError, RateTooLowError +from libecalc.process.process_pipeline.process_error import ( + OutletFluidNotAchievableError, + RateTooHighError, + RateTooLowError, +) from libecalc.process.process_solver.boundary import Boundary from libecalc.process.process_solver.configuration import SpeedConfiguration -from libecalc.process.process_solver.search_strategies import RootFindingStrategy, SearchStrategy +from libecalc.process.process_solver.search_strategies import ( + ACCEPT_AND_GO_HIGHER, + ACCEPT_AND_GO_LOWER, + REJECT_AND_GO_HIGHER, + REJECT_AND_GO_LOWER, + BinarySearchResult, + RootFindingStrategy, + SearchStrategy, +) from libecalc.process.process_solver.solver import ( + OutletFluidNotAchievableEvent, OutsideCapacityEvent, Solution, Solver, @@ -61,6 +74,12 @@ def get_outlet_stream(speed: float) -> FluidStream: source_id=e.process_unit_id, ), ) + except OutletFluidNotAchievableError: + result = self._search_highest_eos_feasible_speed(get_outlet_stream) + if isinstance(result, Solution): + return result + max_speed_configuration = result + maximum_speed_outlet_stream = get_outlet_stream(speed=result.speed) if maximum_speed_outlet_stream.pressure_bara < self._target_pressure: return Solution( @@ -95,6 +114,12 @@ def bool_speed_func(x: float): ) minimum_speed_configuration = SpeedConfiguration(speed=minimum_speed_within_capacity) minimum_speed_outlet_stream = func(minimum_speed_configuration) + except OutletFluidNotAchievableError: + result = self._search_lowest_eos_feasible_speed(get_outlet_stream) + if isinstance(result, Solution): + return result + minimum_speed_configuration = result + minimum_speed_outlet_stream = get_outlet_stream(speed=result.speed) if minimum_speed_outlet_stream.pressure_bara > self._target_pressure: # Solution 2, target pressure is too low @@ -121,8 +146,85 @@ def root_speed_func(x: float) -> float: out = get_outlet_stream(speed=x) return out.pressure_bara - self._target_pressure - speed = self._root_finding_strategy.find_root( - boundary=Boundary(min=minimum_speed_configuration.speed, max=self._boundary.max), - func=root_speed_func, - ) + try: + speed = self._root_finding_strategy.find_root( + boundary=Boundary(min=minimum_speed_configuration.speed, max=max_speed_configuration.speed), + func=root_speed_func, + ) + except OutletFluidNotAchievableError as e: + return self._fluid_not_achievable_solution( + e, SpeedConfiguration(speed=e.unachievable_operating_point.speed) + ) + return Solution(success=True, configuration=SpeedConfiguration(speed=speed)) + + def _fluid_not_achievable_solution( + self, + e: OutletFluidNotAchievableError, + configuration: SpeedConfiguration, + ) -> Solution[SpeedConfiguration]: + logger.warning(f"Outlet fluid not achievable at speed {configuration.speed}", exc_info=e) + return Solution( + success=False, + configuration=configuration, + failure_event=OutletFluidNotAchievableEvent( + status=SolverFailureStatus.OUTLET_FLUID_NOT_ACHIEVABLE, + source_id=e.process_unit_id, + unachievable_operating_point=e.unachievable_operating_point, + ), + ) + + def _search_highest_eos_feasible_speed( + self, + get_outlet_stream: Callable[[float], FluidStream], + ) -> SpeedConfiguration | Solution[SpeedConfiguration]: + """Binary-search downward for the highest speed where the EOS produces a valid outlet. + + Pre-checks that boundary.min is feasible; if not, returns a failed Solution + (the search would otherwise exhaust iterations without converging). + """ + logger.debug( + "Outlet fluid not achievable at maximum speed %.1f — searching for highest feasible speed.", + self._boundary.max, + ) + try: + get_outlet_stream(self._boundary.min) + except OutletFluidNotAchievableError as e: + return self._fluid_not_achievable_solution(e, SpeedConfiguration(speed=self._boundary.min)) + + def feasible(x: float) -> BinarySearchResult: + try: + get_outlet_stream(x) + return ACCEPT_AND_GO_HIGHER + except (OutletFluidNotAchievableError, RateTooHighError, RateTooLowError): + return REJECT_AND_GO_LOWER + + speed = self._search_strategy.search(boundary=self._boundary, func=feasible) + return SpeedConfiguration(speed=speed) + + def _search_lowest_eos_feasible_speed( + self, + get_outlet_stream: Callable[[float], FluidStream], + ) -> SpeedConfiguration | Solution[SpeedConfiguration]: + """Binary-search upward for the lowest speed where the EOS produces a valid outlet. + + Pre-checks that boundary.max is feasible; if not, returns a failed Solution. + """ + logger.debug( + "Outlet fluid not achievable at minimum speed %.1f — searching for lowest feasible speed.", + self._boundary.min, + ) + try: + get_outlet_stream(self._boundary.max) + except OutletFluidNotAchievableError as e: + return self._fluid_not_achievable_solution(e, SpeedConfiguration(speed=self._boundary.max)) + + def feasible(x: float) -> BinarySearchResult: + try: + get_outlet_stream(x) + return ACCEPT_AND_GO_LOWER + except (OutletFluidNotAchievableError, RateTooHighError, RateTooLowError): + return REJECT_AND_GO_HIGHER + + speed = self._search_strategy.search(boundary=self._boundary, func=feasible) + return SpeedConfiguration(speed=speed) diff --git a/src/libecalc/process/process_units/compressor.py b/src/libecalc/process/process_units/compressor.py index 8340a11ab7..e2a9be4eb5 100644 --- a/src/libecalc/process/process_units/compressor.py +++ b/src/libecalc/process/process_units/compressor.py @@ -1,5 +1,6 @@ from typing import Final +from libecalc.domain.process.compressor.core.exceptions import CompressorThermodynamicCalculationError from libecalc.domain.process.compressor.core.train.utils.common import ( RECIRCULATION_BOUNDARY_TOLERANCE, calculate_outlet_pressure_and_stream, @@ -8,7 +9,12 @@ from libecalc.domain.process.value_objects.chart.compressor import CompressorChart from libecalc.process.fluid_stream.fluid_service import FluidService from libecalc.process.fluid_stream.fluid_stream import FluidStream -from libecalc.process.process_pipeline.process_error import RateTooHighError, RateTooLowError +from libecalc.process.process_pipeline.process_error import ( + CompressorOperatingPoint, + OutletFluidNotAchievableError, + RateTooHighError, + RateTooLowError, +) from libecalc.process.process_pipeline.process_unit import ProcessUnit, ProcessUnitId from libecalc.process.process_solver.boundary import Boundary @@ -56,12 +62,27 @@ def propagate_stream(self, inlet_stream: FluidStream) -> FluidStream: rate=actual_rate, ) - return calculate_outlet_pressure_and_stream( - polytropic_efficiency=polytropic_efficiency, - polytropic_head_joule_per_kg=polytropic_head, - inlet_stream=inlet_stream, - fluid_service=self._fluid_service, - ) + try: + return calculate_outlet_pressure_and_stream( + polytropic_efficiency=polytropic_efficiency, + polytropic_head_joule_per_kg=polytropic_head, + inlet_stream=inlet_stream, + fluid_service=self._fluid_service, + ) + except CompressorThermodynamicCalculationError as exc: + # The compressor outlet thermodynamics could not produce a usable state + # (invalid Campbell pressure guess, PH flash failure, or invalid PH result). + raise OutletFluidNotAchievableError( + process_unit_id=self._id, + unachievable_operating_point=CompressorOperatingPoint( + inlet_pressure_bara=inlet_stream.pressure_bara, + inlet_temperature_kelvin=inlet_stream.temperature_kelvin, + actual_rate_m3_per_hour=actual_rate, + polytropic_head_joule_per_kg=polytropic_head, + polytropic_efficiency=polytropic_efficiency, + speed=self.speed, + ), + ) from exc @property def compressor_chart(self) -> CompressorChart: diff --git a/tests/libecalc/process/process_solver/solvers/test_speed_solver.py b/tests/libecalc/process/process_solver/solvers/test_speed_solver.py index 9c38508cbb..9137e168ec 100644 --- a/tests/libecalc/process/process_solver/solvers/test_speed_solver.py +++ b/tests/libecalc/process/process_solver/solvers/test_speed_solver.py @@ -1,12 +1,18 @@ +from collections.abc import Callable from typing import Final import pytest from libecalc.process.fluid_stream.fluid_service import FluidService from libecalc.process.fluid_stream.fluid_stream import FluidStream +from libecalc.process.process_pipeline.process_error import ( + CompressorOperatingPoint, + OutletFluidNotAchievableError, +) from libecalc.process.process_pipeline.process_unit import ProcessUnit, ProcessUnitId from libecalc.process.process_solver.boundary import Boundary from libecalc.process.process_solver.process_pipeline_runner import propagate_stream_many +from libecalc.process.process_solver.solver import SolverFailureStatus, TargetNotAchievableEvent from libecalc.process.process_solver.solvers.speed_solver import SpeedConfiguration, SpeedSolver from libecalc.process.shaft import Shaft, VariableSpeedShaft @@ -85,3 +91,227 @@ def speed_func(configuration: SpeedConfiguration): assert solution.configuration.speed == expected_speed outlet_stream = speed_func(solution.configuration) assert outlet_stream.pressure_bara == expected_pressure + + +class FakeProcessUnit(ProcessUnit): + """Test double whose propagate_stream is supplied as a callable.""" + + def __init__( + self, + propagate_stream: Callable[[FluidStream], FluidStream], + process_unit_id: ProcessUnitId | None = None, + ): + self._id: Final[ProcessUnitId] = process_unit_id or ProcessUnit._create_id() + self._propagate_stream = propagate_stream + + def get_id(self) -> ProcessUnitId: + return self._id + + def propagate_stream(self, inlet_stream: FluidStream) -> FluidStream: + return self._propagate_stream(inlet_stream) + + +def _failing_speed_unit( + shaft: Shaft, + fluid_service: FluidService, + *, + fails_at_or_above: float | None = None, + fails_at_or_below: float | None = None, +) -> FakeProcessUnit: + """Build a FakeProcessUnit that mirrors `pressure = inlet + speed` but raises + OutletFluidNotAchievableError when the shaft speed crosses the given threshold.""" + unit_id = ProcessUnit._create_id() + + def _propagate(inlet_stream: FluidStream) -> FluidStream: + speed = shaft.get_speed() + should_fail = (fails_at_or_above is not None and speed >= fails_at_or_above) or ( + fails_at_or_below is not None and speed <= fails_at_or_below + ) + if should_fail: + raise OutletFluidNotAchievableError( + process_unit_id=unit_id, + unachievable_operating_point=CompressorOperatingPoint( + inlet_pressure_bara=inlet_stream.pressure_bara, + inlet_temperature_kelvin=inlet_stream.temperature_kelvin, + actual_rate_m3_per_hour=0.0, + polytropic_head_joule_per_kg=0.0, + polytropic_efficiency=0.0, + speed=speed, + ), + ) + return fluid_service.create_stream_from_standard_rate( + fluid_model=inlet_stream.fluid_model, + pressure_bara=inlet_stream.pressure_bara + speed, + standard_rate_m3_per_day=inlet_stream.standard_rate_sm3_per_day, + temperature_kelvin=inlet_stream.temperature_kelvin, + ) + + return FakeProcessUnit(propagate_stream=_propagate, process_unit_id=unit_id) + + +def test_min_speed_fluid_not_achievable_target_achievable( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """OutletFluidNotAchievableError at min speed: binary search finds higher effective min; target still reachable.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_below=300) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=450, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is True + assert abs(solution.configuration.speed - 350) < 1.0 + + +def test_min_speed_fluid_not_achievable_target_not_achievable( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """OutletFluidNotAchievableError at min speed: effective min too high → target below minimum achievable.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_below=400) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=350, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is False + assert isinstance(solution.failure_event, TargetNotAchievableEvent) + assert solution.failure_event.status == SolverFailureStatus.MINIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_ABOVE_TARGET + assert solution.failure_event.achievable_value > 350 + assert solution.failure_event.target_value == 350 + + +def test_max_speed_fluid_not_achievable_target_achievable( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """OutletFluidNotAchievableError at max speed: binary search finds lower effective max; target still reachable.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_above=500) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=350, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is True + assert abs(solution.configuration.speed - 250) < 1.0 + + +def test_max_speed_fluid_not_achievable_target_not_achievable( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """OutletFluidNotAchievableError at max speed: binary search finds lower effective max; target not reachable.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_above=300) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=500, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is False + assert isinstance(solution.failure_event, TargetNotAchievableEvent) + assert solution.failure_event.status == SolverFailureStatus.MAXIMUM_ACHIEVABLE_DISCHARGE_PRESSURE_BELOW_TARGET + assert solution.failure_event.achievable_value < 500 + assert solution.failure_event.target_value == 500 + + +def test_all_speeds_fluid_not_achievable_from_max( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """All speeds fail EOS when searching from max: pre-check at boundary.min catches it immediately.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_above=0) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=350, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is False + assert solution.failure_event.status == SolverFailureStatus.OUTLET_FLUID_NOT_ACHIEVABLE + + +def test_all_speeds_fluid_not_achievable_from_min( + search_strategy_factory, + root_finding_strategy, + stream_factory, + shaft, + fluid_service, +): + """All speeds fail EOS when searching from min: pre-check at boundary.max catches it immediately.""" + unit = _failing_speed_unit(shaft=shaft, fluid_service=fluid_service, fails_at_or_below=600) + speed_solver = SpeedSolver( + search_strategy=search_strategy_factory(), + root_finding_strategy=root_finding_strategy, + boundary=Boundary(min=200, max=600), + target_pressure=350, + ) + inlet_stream = stream_factory(standard_rate_m3_per_day=1000, pressure_bara=100) + + def speed_func(configuration: SpeedConfiguration) -> FluidStream: + shaft.set_speed(configuration.speed) + return propagate_stream_many(process_units=[unit], inlet_stream=inlet_stream) + + solution = speed_solver.solve(speed_func) + + assert solution.success is False + assert solution.failure_event.status == SolverFailureStatus.OUTLET_FLUID_NOT_ACHIEVABLE