Add IST reduction

src/thermohl/power/__init__.py

@@ -8,12 +8,13 @@
 from .radiative_cooling import RadiativeCoolingBase
 
 from .power_term import PowerTerm
-from .solar_heating import _SRad, SolarHeatingBase
+from .solar_heating import _SRad, SolarHeatingBase, FixedSolarIrradianceSolarHeating
 
 
 __all__ = [
     "RadiativeCoolingBase",
     "PowerTerm",
     "_SRad",
     "SolarHeatingBase",
+    "FixedSolarIrradianceSolarHeating",
 ]

src/thermohl/power/cigre/convective_cooling.py

@@ -11,15 +11,16 @@
 import numpy as np
 
 from thermohl import floatArrayLike
-from thermohl.power import PowerTerm
-from thermohl.power.convective_cooling import compute_wind_attack_angle
+from thermohl.power.convective_cooling import (
+    ConvectiveCoolingBase,
+)
 from thermohl.power.cigre import Air
 
 
 logger = logging.getLogger(__name__)
 
 
-class ConvectiveCooling(PowerTerm):
+class ConvectiveCooling(ConvectiveCoolingBase):
     """Convective cooling term."""
 
     def __init__(
@@ -30,8 +31,8 @@ def __init__(
         wind_speed: floatArrayLike,
         outer_diameter: floatArrayLike,
         roughness_ratio: floatArrayLike,
-        wind_azimuth: floatArrayLike = None,
-        wind_attack_angle: floatArrayLike = None,
+        wind_azimuth: floatArrayLike | None = None,
+        wind_attack_angle: floatArrayLike | None = None,
         g: float = 9.81,
         **kwargs: Any,
     ):
@@ -44,33 +45,24 @@ def __init__(
             cable_azimuth (float | numpy.ndarray): Azimuth (deg).
             ambient_temperature (float | numpy.ndarray): Ambient temperature (°C).
             wind_speed (float | numpy.ndarray): Wind speed (m·s⁻¹).
-            wind_azimuth (float | numpy.ndarray): wind azimuth regarding north (deg).
+            wind_azimuth (float | numpy.ndarray | None): wind azimuth regarding north (deg).
+            wind_attack_angle (float | numpy.ndarray | None): wind attack angle (rad).
             outer_diameter (float | numpy.ndarray): External diameter (m).
             roughness_ratio (float | numpy.ndarray): Cable roughness (—).
             g (float, optional): Gravitational acceleration (m·s⁻²). The default is 9.81.
 
         """
+        self._check_arguments(wind_azimuth, wind_attack_angle)
+        self._set_wind_attack_angle(cable_azimuth, wind_azimuth, wind_attack_angle)
+
         self.altitude = altitude
         self.ambient_temp = ambient_temperature
         self.wind_speed = wind_speed
         self.outer_diameter = outer_diameter
         self.roughness_ratio = roughness_ratio
         self.gravity = g
 
-        if wind_attack_angle is None and wind_azimuth is None:
-            raise ValueError("Must provide either wind_attack_angle or wind_azimuth.")
-        if wind_attack_angle is not None and wind_azimuth is not None:
-            logger.warning(
-                "both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
-            )
-        if wind_attack_angle is not None:
-            self.wind_attack_angle = wind_attack_angle
-        else:
-            self.wind_attack_angle = compute_wind_attack_angle(
-                cable_azimuth, wind_azimuth
-            )
-
-    def _nu_forced(
+    def _nusselt_forced(
         self, film_temperature: floatArrayLike, kinematic_viscosity: floatArrayLike
     ) -> floatArrayLike:
         """
@@ -133,7 +125,7 @@ def _nu_forced(
             B1 * reynolds**n
         )
 
-    def _nu_natural(
+    def _nusselt_natural(
         self,
         film_temperature: floatArrayLike,
         temperature_delta: floatArrayLike,
@@ -192,8 +184,8 @@ def value(self, conductor_temperature: floatArrayLike) -> floatArrayLike:
         # nu[nu < 1.0E-06] = 1.0E-06
         thermal_conductivity = Air.thermal_conductivity(film_temperature)
         # lm[lm < 0.01] = 0.01
-        nusselt_forced = self._nu_forced(film_temperature, kinematic_viscosity)
-        nusselt_natural = self._nu_natural(
+        nusselt_forced = self._nusselt_forced(film_temperature, kinematic_viscosity)
+        nusselt_natural = self._nusselt_natural(
             film_temperature, temperature_delta, kinematic_viscosity
         )
         return (

src/thermohl/power/convective_cooling.py

@@ -46,32 +46,77 @@ def __init__(
         air_density: Callable[[floatArrayLike, floatArrayLike], floatArrayLike],
         dynamic_viscosity: Callable[[floatArrayLike], floatArrayLike],
         thermal_conductivity: Callable[[floatArrayLike], floatArrayLike],
-        wind_azimuth: floatArrayLike = None,
-        wind_attack_angle: floatArrayLike = None,
+        wind_azimuth: floatArrayLike | None = None,
+        wind_attack_angle: floatArrayLike | None = None,
         **kwargs: Any,
     ):
+        self._check_arguments(wind_azimuth, wind_attack_angle)
+        self._set_wind_attack_angle(cable_azimuth, wind_azimuth, wind_attack_angle)
+
         self.altitude = altitude
         self.ambient_temp = ambient_temperature
         self.wind_speed = wind_speed
 
-        if wind_attack_angle is None and wind_azimuth is None:
+        self.outer_diameter = outer_diameter
+
+        self.air_density = air_density
+        self.dynamic_viscosity = dynamic_viscosity
+        self.thermal_conductivity = thermal_conductivity
+
+    @classmethod
+    def _check_arguments(
+        cls,
+        wind_azimuth: floatArrayLike | None,
+        wind_attack_angle: floatArrayLike | None,
+    ) -> None:
+        if (wind_attack_angle is None or np.isnan(wind_attack_angle).any()) and (
+            wind_azimuth is None or np.isnan(wind_azimuth).any()
+        ):
             raise ValueError("Must provide either wind_attack_angle or wind_azimuth.")
-        if wind_attack_angle is not None and wind_azimuth is not None:
+
+        if (
+            wind_attack_angle is not None
+            and not np.isnan(wind_attack_angle).all()
+            and wind_azimuth is not None
+        ):
             logger.warning(
-                "both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
+                "Both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored."
             )
-        if wind_attack_angle is not None:
+
+    def _set_wind_attack_angle(
+        self,
+        cable_azimuth: floatArrayLike,
+        wind_azimuth: floatArrayLike,
+        wind_attack_angle: floatArrayLike,
+    ) -> None:
+        # Compute missing wind attack angles
+        if isinstance(wind_attack_angle, np.ndarray) and wind_attack_angle.ndim > 0:
+            if (
+                not isinstance(cable_azimuth, np.ndarray)
+                or cable_azimuth.shape != wind_attack_angle.shape
+            ):
+                raise ValueError(
+                    "If wind_attack_angle is an array, cable_azimuth must be an array of the same shape."
+                )
+            if (
+                not isinstance(wind_azimuth, np.ndarray)
+                or wind_azimuth.shape != wind_attack_angle.shape
+            ):
+                raise ValueError(
+                    "If wind_attack_angle is an array, wind_azimuth must be an array of the same shape."
+                )
+            mask = np.isnan(wind_attack_angle)
+            if np.any(mask):
+                wind_attack_angle[mask] = compute_wind_attack_angle(
+                    cable_azimuth[mask], wind_azimuth[mask]
+                )
             self.wind_attack_angle = wind_attack_angle
-        else:
+        elif wind_attack_angle is None or np.isnan(wind_attack_angle):
             self.wind_attack_angle = compute_wind_attack_angle(
                 cable_azimuth, wind_azimuth
             )
-
-        self.outer_diameter = outer_diameter
-
-        self.air_density = air_density
-        self.dynamic_viscosity = dynamic_viscosity
-        self.thermal_conductivity = thermal_conductivity
+        else:
+            self.wind_attack_angle = wind_attack_angle
 
     def _value_forced(
         self,

src/thermohl/power/solar_heating.py

@@ -147,3 +147,21 @@ def derivative(self, conductor_temperature: floatArrayLike) -> floatArrayLike:
         :return: Derivative of solar heating.
         """
         return np.zeros_like(conductor_temperature)
+
+
+class FixedSolarIrradianceSolarHeating(SolarHeatingBase):
+    """Solar heating term with fixed solar irradiance.
+
+    This class only computes the solar heating (power) and its
+    derivative based on the provided solar irradiance.
+    """
+
+    def __init__(
+        self,
+        outer_diameter: floatArrayLike,
+        solar_absorptivity: floatArrayLike,
+        solar_irradiance: floatArrayLike,
+    ):
+        self.outer_diameter = outer_diameter
+        self.solar_absorptivity = solar_absorptivity
+        self.solar_irradiance = solar_irradiance

src/thermohl/solver/parameters.py

@@ -51,6 +51,7 @@ def _set_default_values(self) -> None:
         self.precipitation_rate = 0.0  # rain precipitation rate (m.s**-1)
         self.wind_speed = 0.0  # wind speed (m.s**-1)
         self.wind_azimuth = 90.0  # wind_azimuth (deg, regarding north)
+        self.wind_attack_angle = np.nan  # wind attack angle (rad)
         self.nebulosity = np.nan  # nebulosity (1)
         self.albedo = 0.15  # albedo (1)
         # coefficient for air pollution from 0 (clean) to 1 (polluted)

src/thermohl/solver/slv1t.py

@@ -5,16 +5,22 @@
 # file, You can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 
+import logging
 import numbers
 from typing import Optional
 
 import numpy as np
 
 from thermohl import floatArrayLike, floatArray
+from thermohl.power import FixedSolarIrradianceSolarHeating
 from thermohl.solver.solver import Solver as Solver_, get_time_changing_parameters
 from thermohl.solver.parameters import DEFAULT_PARAMETERS as default
 from thermohl.solver.entities import PowerType, VariableType
 from thermohl.utils import bisect_v
+from thermohl.utils import quasi_newton
+
+
+logger = logging.getLogger(__name__)
 
 
 class Solver1T(Solver_):
@@ -220,6 +226,112 @@ def fun(i: floatArray) -> floatArrayLike:
             return_power,
         )
 
-        result = self._add_input_data_to_result(result)
+        return self._add_input_data_to_result(result)
 
-        return result
+    def _set_default_reduced_intensity_args(
+        self,
+        ambient_temperature: Optional[floatArrayLike],
+        wind_speed: Optional[floatArrayLike],
+        solar_irradiance: Optional[floatArrayLike],
+    ):
+        if ambient_temperature is None:
+            logger.warning(
+                "ambient_temperature is not set. Using default value of 30 °C."
+            )
+            ambient_temperature = 30.0
+        self.args.ambient_temperature = ambient_temperature
+
+        if wind_speed is None:
+            logger.warning("wind_speed is not set. Using default value of 0.6 m/s.")
+            wind_speed = 0.6
+        self.args.wind_speed = wind_speed
+
+        if solar_irradiance is None:
+            logger.warning(
+                "solar_irradiance is not set. Using default value of 600 W/m²."
+            )
+            solar_irradiance = 600.0
+        return solar_irradiance
+
+    def reduced_intensity(
+        self,
+        measured_temperature_difference: floatArrayLike,
+        measured_intensity: floatArrayLike,
+        ambient_temperature: Optional[floatArrayLike] = None,
+        wind_speed: Optional[floatArrayLike] = None,
+        solar_irradiance: Optional[floatArrayLike] = None,
+        max_conductor_temperature: Optional[floatArrayLike] = None,
+    ) -> floatArrayLike:
+        """
+        Compute the reduced intensity limit for a given measured temperature difference
+        betwwen the sound cable and a hotspot on the junction between a cable
+        and a faulty sleeve.
+
+        Args:
+            measured_temperature_difference (float | np.ndarray): The measured temperature difference between the cable surface and the sleeve.
+            measured_intensity (float | np.ndarray): The measuredintensity at which the temperature difference was measured.
+            ambient_temperature (Optional[float | np.ndarray]): The ambient temperature. Default is 30.
+            wind_speed (Optional[float | np.ndarray]): The wind speed. Default is 0.6.
+            solar_irradiance (Optional[float | np.ndarray]): The measured solar irradiance. Default is 600.
+            max_conductor_temperature (Optional[float | np.ndarray]): The maximum conductor temperature. Default is 100.
+        """
+        # Save elements that will be modified to do the computations
+        # so as to be able to restore afterwards
+        saved_solver_transit = self.args.transit
+        saved_solver_ambient_temperature = self.args.ambient_temperature
+        saved_solver_wind_speed = self.args.wind_speed
+        saved_solver_wind_attack_angle = self.args.wind_attack_angle
+
+        saved_solar_heating = self.solar_heating
+
+        # Set args default values for reduced intensity computation.
+        # These differ from those used for the other computations.
+        solar_irradiance = self._set_default_reduced_intensity_args(
+            ambient_temperature, wind_speed, solar_irradiance
+        )
+
+        # Set default value for max_conductor_temperature if not provided.
+        if max_conductor_temperature is None:
+            max_conductor_temperature = np.full_like(measured_intensity, 100.0)
+
+        self.args.wind_attack_angle = 90.0
+        self.convective_cooling.__init__(**self.args.__dict__)
+
+        self.solar_heating = FixedSolarIrradianceSolarHeating(
+            outer_diameter=self.args.outer_diameter,
+            solar_absorptivity=self.args.solar_absorptivity,
+            solar_irradiance=solar_irradiance,
+        )
+
+        def conductor_temperature(transit):
+            self.args.transit = transit
+            self.joule_heating.__init__(**self.args.__dict__)
+            return self.steady_temperature()[VariableType.TEMPERATURE][0]
+
+        def temperature_difference(transit):
+            return measured_temperature_difference * (
+                (transit / measured_intensity) ** 2
+            )
+
+        def sleeve_temperature(transit):
+            return conductor_temperature(transit) + temperature_difference(transit)
+
+        def f(transit):
+            return sleeve_temperature(transit) - max_conductor_temperature
+
+        x0 = np.full_like(measured_intensity, 100.0)
+
+        reduced_intensity = quasi_newton(f, x0=x0)
+
+        # Restore previous elements
+        self.args.transit = saved_solver_transit
+        # Update joule heating with restored transit
+        self.joule_heating.__init__(**self.args.__dict__)
+        self.args.ambient_temperature = saved_solver_ambient_temperature
+        self.args.wind_speed = saved_solver_wind_speed
+        self.args.wind_attack_angle = saved_solver_wind_attack_angle
+        # Update convective cooling with restored wind_attack_angle
+        self.convective_cooling.__init__(**self.args.__dict__)
+        self.solar_heating = saved_solar_heating
+
+        return reduced_intensity