From b6951ff105a4dff11cc927a20be12177dadcb048 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Tue, 3 Mar 2026 12:33:38 +0100 Subject: [PATCH 01/16] feat(#97): Add IST reduction feature Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 88 ++++++++++++++++++++++++++++++++++ test/unit/solver/test_slv1t.py | 19 ++++++++ 2 files changed, 107 insertions(+) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index f6868942..facc5e93 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -237,3 +237,91 @@ def fun(i: floatArray) -> floatArrayLike: ) return df + + def _set_default_reduced_intensity_args( + self, + ambient_temperature: Optional[floatArrayLike], + wind_speed: Optional[floatArrayLike], + measured_solar_irradiance: Optional[floatArrayLike], + ): + if ambient_temperature is None: + print( + "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: + print("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 measured_solar_irradiance is None: + print( + "WARNING: measured_solar_irradiance is not set. Using default value of 600 W/m²." + ) + measured_solar_irradiance = 600.0 + self.args.measured_solar_irradiance = measured_solar_irradiance + + def reduced_intensity( + self, + delta_T_measured: floatArrayLike, + measured_intensity: floatArrayLike, + ambient_temperature: Optional[floatArrayLike] = None, + wind_speed: Optional[floatArrayLike] = None, + measured_solar_irradiance: Optional[floatArrayLike] = None, + T_limit: Optional[floatArrayLike] = np.array([100]), + ): + """ + 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: + delta_T_measured (float): The measured temperature difference between the cable surface and the sleeve. + measured_intensity (float): The measuredintensity at which the temperature difference was measured. + ambient_temperature (Optional[float]): The ambient temperature. Default is 30. + wind_speed (Optional[float]): The wind speed. Default is 0.6. + measured_solar_irradiance (Optional[float]): The measured solar irradiance. Default is 600. + T_limit (Optional[float]): The maximum conductor temperature. Default is 100. + """ + # Save args that will be modified so as to be able to restore them at the end of the computation + solver_transit = self.args.transit + solver_ambient_temperature = self.args.ambient_temperature + solver_wind_speed = self.args.wind_speed + solver_measured_solar_irradiance = self.args.measured_solar_irradiance + + # Set default values for reduced intensity computation. + # These differ from those used for the other computations. + self._set_default_reduced_intensity_args( + ambient_temperature, wind_speed, measured_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 delta_T(transit): + return delta_T_measured * ((transit / measured_intensity) ** 2) + + def sleeve_temperature(transit): + return conductor_temperature(transit) + delta_T(transit) + + def f(transit): + return sleeve_temperature(transit) - T_limit + + imax = 4500.0 # Used as upper bound for the bisection method. + # It is a very high value that should not be reached in practice. + + reduced_intensity, _ = bisect_v( + f, DP.imin, imax, output_shape=(self.args.max_len(),) + ) + + # Restore previous args + self.args.transit = solver_transit + self.args.ambient_temperature = solver_ambient_temperature + self.args.wind_speed = solver_wind_speed + self.args.measured_solar_irradiance = solver_measured_solar_irradiance + + return reduced_intensity diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 866a94af..3f764b54 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -203,3 +203,22 @@ def test_steady_intensity_custom_params(solver): assert isinstance(result, pd.DataFrame) assert VariableType.TRANSIT in result.columns + + +def test_reduced_intensity(solver): + solver_args = solver.args.__dict__.copy() + + result = solver.reduced_intensity( + delta_T_measured=10.0, + measured_intensity=360.0, + T_limit=100.0, + ) + + assert isinstance(result, np.ndarray) + assert not np.isnan(result[0]) + + # Check that solver args have not been changed + assert np.isnan(solver.args.measured_solar_irradiance) + assert solver.args.ambient_temperature == solver_args["ambient_temperature"] + assert solver.args.wind_speed == solver_args["wind_speed"] + assert solver.args.transit == solver_args["transit"] From afd5121f85681d28e6f37420b14e300446e2eff5 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Tue, 3 Mar 2026 09:32:46 +0100 Subject: [PATCH 02/16] docs: Add a notebook with an example Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- .../docs/examples/ex_ist_reduction.ipynb | 102 ++++++++++++++++++ thermohl-docs/mkdocs.yml | 9 +- 2 files changed, 107 insertions(+), 4 deletions(-) create mode 100644 thermohl-docs/docs/examples/ex_ist_reduction.ipynb diff --git a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb new file mode 100644 index 00000000..458dfeea --- /dev/null +++ b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb @@ -0,0 +1,102 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 7, + "id": "8c937b14", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from thermohl import solver\n", + "from thermohl.solver.enums.heat_equation_type import HeatEquationType\n", + "from thermohl.solver.enums.variable_type import VariableType\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "81557e4b", + "metadata": {}, + "outputs": [], + "source": [ + "args = {\n", + " \"max_len\": lambda: 1,\n", + " VariableType.TRANSIT: np.array([0]),\n", + " \"ambient_temperature\": np.array([25]),\n", + " \"wind_speed\": np.array([0]),\n", + " \"wind_azimuth\": np.array([0]),\n", + " \"ambient_pressure\": np.array([101325]),\n", + " \"relative_humidity\": np.array([50]),\n", + " \"precipitation_rate\": np.array([0]),\n", + " \"linear_mass\": 1.0,\n", + " \"heat_capacity\": 1.0,\n", + " \"month\": 1,\n", + " \"day\": 1,\n", + " \"hour\": 0,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "0cf8da43", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING: ambient_temperature is not set. Using default value of 30 °C.\n", + "WARNING: wind_speed is not set. Using default value of 0.6 m/s.\n", + "WARNING: measured_solar_irradiance is not set. Using default value of 600 W/m².\n" + ] + }, + { + "data": { + "text/plain": [ + "array([830.11542031])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "slvr = solver.ieee(\n", + " args,\n", + " heat_equation=HeatEquationType.WITH_ONE_TEMPERATURE,\n", + ")\n", + "\n", + "slvr.reduced_intensity(\n", + " delta_T_measured=10.0,\n", + " measured_intensity=360.0,\n", + " T_limit=100.0,\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "thermohl (3.13.12)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/thermohl-docs/mkdocs.yml b/thermohl-docs/mkdocs.yml index 7e4159c9..1a4dd1bd 100644 --- a/thermohl-docs/mkdocs.yml +++ b/thermohl-docs/mkdocs.yml @@ -36,10 +36,11 @@ nav: - Home: index.md - Getting Started: getting-started.md - Examples: - - Example of steady temperature with 1t solver and IEEE model: examples/ex_steady_1t.ipynb - - Example of steady temperature with 3t solver and CIGRE model: examples/ex_steady_3t.ipynb - - Example of transient temperature with 3t solver and OLLA model: examples/ex_transient_3t.ipynb - - Example of steady temperature with 1t solver, all models and all parameters: examples/ex_all_params.ipynb + - Steady temperature with 1t solver and IEEE model: examples/ex_steady_1t.ipynb + - Steady temperature with 3t solver and CIGRE model: examples/ex_steady_3t.ipynb + - Transient temperature with 3t solver and OLLA model: examples/ex_transient_3t.ipynb + - Steady temperature with 1t solver, all models and all parameters: examples/ex_all_params.ipynb + - IST reduction: examples/ex_ist_reduction.ipynb - Parameters: api-reference/parameters.md - User Guide: user-guide.md - Developer Guide: From a5dbd8c3f212bd1d09129de9ef819549715394af Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Tue, 3 Mar 2026 12:19:50 +0100 Subject: [PATCH 03/16] docs: Document IST reduction Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- thermohl-docs/docs/user-guide.md | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/thermohl-docs/docs/user-guide.md b/thermohl-docs/docs/user-guide.md index a6cf2fa2..90965000 100644 --- a/thermohl-docs/docs/user-guide.md +++ b/thermohl-docs/docs/user-guide.md @@ -153,3 +153,27 @@ specific model with three temperatures for the conductor : * the surface temperature; * the average temperature; * the core temperature. + +## IST reduction (intensity limit reduction) + +For most computations, cable sleeves are ignored. +However, if a sleeve is faulty, it can overheat. ThermoHL enables the user to +compute a reduced intensity limit (called reduced IST) so that the interface +between the sleeve and the cable doesn't exceed a given maximum temperature. +Input data are +- the mesaured temperature difference between the hotspot (at the interface between sleeve and cable) +and the cable far away from the sleeve +- the measured transit. + +This computation uses the single temperature model. + +Please note that following default values differ from the defaults for other +computations: + +| Parameter | Default value | +|----|----| +| ambient_temperature | 30.0 °C | +| wind_speed | 0.6 m/s | +| measured_solar_irradiance | 600.0 W/m² | + +TODO: precision on wind_speed and wind_angle if needed From 8f677b511206a661b6f6801f4e68566902ca095e Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Thu, 5 Mar 2026 11:16:32 +0100 Subject: [PATCH 04/16] Use secant (quasi-newton) method Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 14 ++-- src/thermohl/utils.py | 116 ++++++++++++++++++++++++++++++++ test/unit/solver/test_slv1t.py | 117 +++++++++++++++++++++++++++++---- 3 files changed, 229 insertions(+), 18 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index facc5e93..6f1cd2d8 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -16,6 +16,7 @@ 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 bisect_v, quasi_newton class Solver1T(Solver_): @@ -270,7 +271,7 @@ def reduced_intensity( ambient_temperature: Optional[floatArrayLike] = None, wind_speed: Optional[floatArrayLike] = None, measured_solar_irradiance: Optional[floatArrayLike] = None, - T_limit: Optional[floatArrayLike] = np.array([100]), + T_limit: Optional[floatArrayLike] = None, ): """ Compute the reduced intensity limit for a given measured temperature difference @@ -297,6 +298,10 @@ def reduced_intensity( ambient_temperature, wind_speed, measured_solar_irradiance ) + # Set default value for T_limit if not provided. + if T_limit is None: + T_limit = np.ones_like(measured_intensity) * 100.0 + def conductor_temperature(transit): self.args.transit = transit self.joule_heating.__init__(**self.args.__dict__) @@ -311,12 +316,9 @@ def sleeve_temperature(transit): def f(transit): return sleeve_temperature(transit) - T_limit - imax = 4500.0 # Used as upper bound for the bisection method. - # It is a very high value that should not be reached in practice. + x0 = np.ones_like(measured_intensity) * 100 - reduced_intensity, _ = bisect_v( - f, DP.imin, imax, output_shape=(self.args.max_len(),) - ) + reduced_intensity = quasi_newton(f, x0=x0) # Restore previous args self.args.transit = solver_transit diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index 8f388ade..0ef7dd0c 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -12,6 +12,9 @@ import os from functools import wraps from importlib.util import find_spec +import operator +import warnings + import numpy as np import pandas as pd @@ -186,6 +189,119 @@ def bisect_v( return midpoint, abs_error +def _array_quasi_newton(func, x0, tol, maxiter): + """ + A vectorized version of secant method for arrays. + + Do not use this method directly. This method is called from `quasi_newton` + when ``np.size(x0) > 1`` is ``True``. For docstring, see `quasi_newton`. + + Heavily inspired by the implementation of the SciPy library. + """ + # Explicitly copy `x0` as `p` will be modified inplace, but the + # user's array should not be altered. + p = np.array(x0, copy=True) + + failures = np.ones_like(p, dtype=bool) + nz_der = np.ones_like(failures) + + dx = np.finfo(float).eps ** 0.33 + p1 = p * (1 + dx) + np.where(p >= 0, dx, -dx) + q0 = np.asarray(func(p)) + q1 = np.asarray(func(p1)) + active = np.ones_like(p, dtype=bool) + for _ in range(maxiter): + nz_der = q1 != q0 + # stop iterating if all derivatives are zero + if not nz_der.any(): + p = (p1 + p) / 2.0 + break + # Secant Step + dp = (q1 * (p1 - p))[nz_der] / (q1 - q0)[nz_der] + # only update nonzero derivatives + p = np.asarray(p, dtype=np.result_type(p, p1, dp, np.float64)) + p[nz_der] = p1[nz_der] - dp + active_zero_der = ~nz_der & active + p[active_zero_der] = (p1 + p)[active_zero_der] / 2.0 + active &= nz_der # don't assign zero derivatives again + failures[nz_der] = np.abs(dp) >= tol # not yet converged + # stop iterating if there aren't any failures, not incl zero der + if not failures[nz_der].any(): + break + p1, p = p, p1 + q0 = q1 + q1 = np.asarray(func(p1)) + + zero_der = ~nz_der & failures # don't include converged with zero-ders + if zero_der.any(): + nonzero_dp = p1 != p + # non-zero dp, but infinite newton step + zero_der_nz_dp = zero_der & nonzero_dp + if zero_der_nz_dp.any(): + rms = np.sqrt(sum((p1[zero_der_nz_dp] - p[zero_der_nz_dp]) ** 2)) + warnings.warn(f"RMS of {rms:g} reached", RuntimeWarning, stacklevel=3) + elif failures.any(): + all_or_some = "all" if failures.all() else "some" + msg = f"{all_or_some:s} failed to converge after {maxiter:d} iterations" + if failures.all(): + raise RuntimeError(msg) + warnings.warn(msg, RuntimeWarning, stacklevel=3) + + return p + + +def quasi_newton(func, x0, tol=1.48e-8, maxiter=50, rtol=0.0): # TODO: add type hints + """Find the zero of a function using the quasi-Newton (secant) method. + + Heavily inspired by the implementation of optimize.newton in the SciPy library. + """ + if tol <= 0: + raise ValueError(f"tol too small ({tol:g} <= 0)") + maxiter = operator.index(maxiter) + if maxiter < 1: + raise ValueError("maxiter must be greater than 0") + if np.size(x0) > 1: + return _array_quasi_newton(func, x0, tol, maxiter) + + # Convert to float (don't use float(x0); this works also for complex x0) + # Use np.asarray because we want x0 to be a numpy object, not a Python + # object. e.g. np.complex(1+1j) > 0 is possible, but (1 + 1j) > 0 raises + # a TypeError + x0 = np.asarray(x0)[()] * 1.0 + p0 = x0 + + eps = 1e-4 + p1 = x0 * (1 + eps) + p1 += eps if p1 >= 0 else -eps + q0 = func(p0) + q1 = func(p1) + if abs(q1) < abs(q0): + p0, p1, q0, q1 = p1, p0, q1, q0 + for itr in range(maxiter): + if q1 == q0: + if p1 != p0: + msg = f"Tolerance of {p1 - p0} reached." + msg += ( + f" Failed to converge after {itr + 1} iterations," + f" value is {p1}." + ) + raise RuntimeError(msg) + return (p1 + p0) / 2.0 + else: + if abs(q1) > abs(q0): + p = (-q0 / q1 * p1 + p0) / (1 - q0 / q1) + else: + p = (-q1 / q0 * p0 + p1) / (1 - q1 / q0) + if np.isclose(p, p1, rtol=rtol, atol=tol): + return p + p0, q0 = p1, q1 + p1 = p + q1 = func(p1) + + msg = f"Failed to converge after {itr + 1} iterations, value is {p}." + raise RuntimeError(msg) + + # In agreement with Eurobios, this function has been retrieved from the pyntb library, # in order to remove the external dependency on this library. # In this library, this function was initially developed under the name qnewt2d_v diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 3f764b54..a36b42aa 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -16,7 +16,7 @@ @pytest.fixture -def solver(): +def solver_args(): args = { "max_len": lambda: 1, "transit": np.array([0]), @@ -26,8 +26,6 @@ def solver(): "ambient_pressure": np.array([101325]), "relative_humidity": np.array([50]), "precipitation_rate": np.array([0]), - "linear_mass": 1.0, - "heat_capacity": 1.0, "datetime_utc": datetime(2025, 1, 1, 0, tzinfo=timezone.utc), "latitude": np.array([48.0]), "longitude": np.array([2.3]), @@ -36,6 +34,8 @@ def solver(): "solar_absorptivity": np.array([0.5]), } cable = { + "linear_mass": 1.0, + "heat_capacity": 1.0, "outer_diameter": np.array([3.186e-02]), "core_diameter": np.array([0.000]), "outer_area": np.array([6.004e-4]), @@ -48,14 +48,26 @@ def solver(): "emissivity": np.array([0.8]), } args.update(cable) + return args + + +@pytest.fixture +def solver(solver_args): + return create_solver(solver_args) + +def create_solver(solver_args): joule = power.rte.joule_heating.JouleHeating solar = power.rte.solar_heating.SolarHeating convective = power.rte.convective_cooling.ConvectiveCooling radiative = power.rte.radiative_cooling.RadiativeCooling solver = Solver1T( - dic=args, joule=joule, solar=solar, convective=convective, radiative=radiative + dic=solver_args, + joule=joule, + solar=solar, + convective=convective, + radiative=radiative, ) return solver @@ -205,20 +217,101 @@ def test_steady_intensity_custom_params(solver): assert VariableType.TRANSIT in result.columns -def test_reduced_intensity(solver): - solver_args = solver.args.__dict__.copy() +def test_reduced_intensity_scalar(solver): + args = { + "max_len": lambda: 1, + VariableType.TRANSIT.value: 0, + "ambient_temperature": 25, + "wind_speed": 0, + "wind_azimuth": 0, + "ambient_pressure": 101325, + "relative_humidity": 50, + "precipitation_rate": 0, + "month": 1, + "day": 1, + "hour": 0, + "latitude": 48.0, + "longitude": 2.3, + "altitude": 50.0, + "cable_azimuth": 0.0, + "solar_absorptivity": 0.5, + } + cable = { + "linear_mass": 1.0, + "heat_capacity": 1.0, + "outer_diameter": 3.186e-02, + "core_diameter": 0.000, + "outer_area": 6.004e-4, + "core_area": 0.000, + "magnetic_coeff": 1.000, + "magnetic_coeff_per_a": 0.000, + "temperature_coeff_linear": 3.600e-3, + "temperature_coeff_quadratic": 8.000e-7, + "linear_resistance_dc_20c": 5.540e-5, + "emissivity": 0.8, + } + args.update(cable) + solver = create_solver(args) result = solver.reduced_intensity( delta_T_measured=10.0, measured_intensity=360.0, - T_limit=100.0, ) - assert isinstance(result, np.ndarray) - assert not np.isnan(result[0]) + assert isinstance(result, np.float64) + assert not np.isnan(result) # Check that solver args have not been changed assert np.isnan(solver.args.measured_solar_irradiance) - assert solver.args.ambient_temperature == solver_args["ambient_temperature"] - assert solver.args.wind_speed == solver_args["wind_speed"] - assert solver.args.transit == solver_args["transit"] + assert solver.args.ambient_temperature == args["ambient_temperature"] + assert solver.args.wind_speed == args["wind_speed"] + assert solver.args.transit == args["transit"] + + +def test_reduced_intensity_array(): + args = { + "max_len": lambda: 2, + VariableType.TRANSIT.value: np.array([0, 0]), + "ambient_temperature": np.array([25, 25]), + "wind_speed": np.array([0, 0]), + "wind_azimuth": np.array([0, 0]), + "ambient_pressure": np.array([101325, 101325]), + "relative_humidity": np.array([50, 50]), + "precipitation_rate": np.array([0, 0]), + "month": np.array([1, 1]), + "day": np.array([1, 1]), + "hour": np.array([0, 0]), + "latitude": np.array([48.0, 48.0]), + "longitude": np.array([2.3, 2.3]), + "altitude": np.array([50.0, 50.0]), + "cable_azimuth": np.array([0.0, 0.0]), + "solar_absorptivity": np.array([0.5, 0.5]), + } + cable = { + "linear_mass": np.array([1.0, 1.0]), + "heat_capacity": np.array([1.0, 1.0]), + "outer_diameter": np.array([3.186e-02, 3.186e-02]), + "core_diameter": np.array([0.000, 0.000]), + "outer_area": np.array([6.004e-4, 6.004e-4]), + "core_area": np.array([0.000, 0.000]), + "magnetic_coeff": np.array([1.000, 1.000]), + "magnetic_coeff_per_a": np.array([0.000, 0.000]), + "temperature_coeff_linear": np.array([3.600e-3, 3.600e-3]), + "temperature_coeff_quadratic": np.array([8.000e-7, 8.000e-7]), + "linear_resistance_dc_20c": np.array([5.540e-5, 5.540e-5]), + "emissivity": np.array([0.8, 0.8]), + } + args.update(cable) + + solver = create_solver(args) + + result = solver.reduced_intensity( + delta_T_measured=np.array([10.0, 10.0]), + measured_intensity=np.array([360.0, 360.0]), + T_limit=np.array([100.0, 100.0]), + ) + + assert isinstance(result, np.ndarray) + assert len(result) == 2 + assert not np.isnan(result[0]) + assert not np.isnan(result[1]) From bdfd1816dae89da271d9904aa1e0f62169692263 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Mon, 9 Mar 2026 11:10:34 +0100 Subject: [PATCH 05/16] Fix some sonarqube errors Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 24 ++++++++------- src/thermohl/utils.py | 30 +++++++++++++------ test/unit/solver/test_slv1t.py | 8 ++--- .../docs/examples/ex_ist_reduction.ipynb | 14 ++++----- 4 files changed, 44 insertions(+), 32 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 6f1cd2d8..eca9d131 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -266,12 +266,12 @@ def _set_default_reduced_intensity_args( def reduced_intensity( self, - delta_T_measured: floatArrayLike, + measured_temperature_difference: floatArrayLike, measured_intensity: floatArrayLike, ambient_temperature: Optional[floatArrayLike] = None, wind_speed: Optional[floatArrayLike] = None, measured_solar_irradiance: Optional[floatArrayLike] = None, - T_limit: Optional[floatArrayLike] = None, + max_conductor_temperature: Optional[floatArrayLike] = None, ): """ Compute the reduced intensity limit for a given measured temperature difference @@ -279,12 +279,12 @@ def reduced_intensity( and a faulty sleeve. Args: - delta_T_measured (float): The measured temperature difference between the cable surface and the sleeve. + measured_temperature_difference (float): The measured temperature difference between the cable surface and the sleeve. measured_intensity (float): The measuredintensity at which the temperature difference was measured. ambient_temperature (Optional[float]): The ambient temperature. Default is 30. wind_speed (Optional[float]): The wind speed. Default is 0.6. measured_solar_irradiance (Optional[float]): The measured solar irradiance. Default is 600. - T_limit (Optional[float]): The maximum conductor temperature. Default is 100. + max_conductor_temperature (Optional[float]): The maximum conductor temperature. Default is 100. """ # Save args that will be modified so as to be able to restore them at the end of the computation solver_transit = self.args.transit @@ -298,23 +298,25 @@ def reduced_intensity( ambient_temperature, wind_speed, measured_solar_irradiance ) - # Set default value for T_limit if not provided. - if T_limit is None: - T_limit = np.ones_like(measured_intensity) * 100.0 + # Set default value for max_conductor_temperature if not provided. + if max_conductor_temperature is None: + max_conductor_temperature = np.ones_like(measured_intensity) * 100.0 def conductor_temperature(transit): self.args.transit = transit self.joule_heating.__init__(**self.args.__dict__) return self.steady_temperature()[VariableType.TEMPERATURE][0] - def delta_T(transit): - return delta_T_measured * ((transit / measured_intensity) ** 2) + def temperature_difference(transit): + return measured_temperature_difference * ( + (transit / measured_intensity) ** 2 + ) def sleeve_temperature(transit): - return conductor_temperature(transit) + delta_T(transit) + return conductor_temperature(transit) + temperature_difference(transit) def f(transit): - return sleeve_temperature(transit) - T_limit + return sleeve_temperature(transit) - max_conductor_temperature x0 = np.ones_like(measured_intensity) * 100 diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index 0ef7dd0c..d96e11f2 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -12,14 +12,15 @@ import os from functools import wraps from importlib.util import find_spec -import operator import warnings - +from typing import Callable import numpy as np import pandas as pd import yaml +from thermohl import floatArrayLike + def _dict_completion( params: dict, @@ -189,7 +190,9 @@ def bisect_v( return midpoint, abs_error -def _array_quasi_newton(func, x0, tol, maxiter): +def _array_quasi_newton( + func: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, tol: float, maxiter: int +) -> floatArrayLike: """ A vectorized version of secant method for arrays. @@ -250,16 +253,25 @@ def _array_quasi_newton(func, x0, tol, maxiter): return p -def quasi_newton(func, x0, tol=1.48e-8, maxiter=50, rtol=0.0): # TODO: add type hints - """Find the zero of a function using the quasi-Newton (secant) method. - - Heavily inspired by the implementation of optimize.newton in the SciPy library. - """ +def _check_quasi_newton_arguments(tol: float, maxiter: int) -> None: if tol <= 0: raise ValueError(f"tol too small ({tol:g} <= 0)") - maxiter = operator.index(maxiter) if maxiter < 1: raise ValueError("maxiter must be greater than 0") + + +def quasi_newton( + func: Callable[[floatArrayLike], floatArrayLike], + x0: floatArrayLike, + tol: float = 1.48e-8, + maxiter: int = 50, + rtol: float = 0.0, +) -> floatArrayLike: + """Find the zero of a function using the quasi-Newton (secant) method. + + Heavily inspired by the implementation of optimize.newton in the SciPy library. + """ + _check_quasi_newton_arguments(tol, maxiter) if np.size(x0) > 1: return _array_quasi_newton(func, x0, tol, maxiter) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index a36b42aa..540e31cd 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -217,7 +217,7 @@ def test_steady_intensity_custom_params(solver): assert VariableType.TRANSIT in result.columns -def test_reduced_intensity_scalar(solver): +def test_reduced_intensity_scalar(): args = { "max_len": lambda: 1, VariableType.TRANSIT.value: 0, @@ -254,12 +254,11 @@ def test_reduced_intensity_scalar(solver): solver = create_solver(args) result = solver.reduced_intensity( - delta_T_measured=10.0, + measured_temperature_difference=10.0, measured_intensity=360.0, ) assert isinstance(result, np.float64) - assert not np.isnan(result) # Check that solver args have not been changed assert np.isnan(solver.args.measured_solar_irradiance) @@ -306,9 +305,8 @@ def test_reduced_intensity_array(): solver = create_solver(args) result = solver.reduced_intensity( - delta_T_measured=np.array([10.0, 10.0]), + measured_temperature_difference=np.array([10.0, 10.0]), measured_intensity=np.array([360.0, 360.0]), - T_limit=np.array([100.0, 100.0]), ) assert isinstance(result, np.ndarray) diff --git a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb index 458dfeea..789a16af 100644 --- a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb +++ b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 7, + "execution_count": 1, "id": "8c937b14", "metadata": {}, "outputs": [], @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 2, "id": "81557e4b", "metadata": {}, "outputs": [], @@ -40,7 +40,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 3, "id": "0cf8da43", "metadata": {}, "outputs": [ @@ -56,10 +56,10 @@ { "data": { "text/plain": [ - "array([830.11542031])" + "np.float64(830.1154202065042)" ] }, - "execution_count": 9, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -71,9 +71,9 @@ ")\n", "\n", "slvr.reduced_intensity(\n", - " delta_T_measured=10.0,\n", + " measured_temperature_difference=10.0,\n", " measured_intensity=360.0,\n", - " T_limit=100.0,\n", + " max_conductor_temperature=100.0,\n", ")" ] } From a85d00f2aeba4c2eae75c05a5f74288e1fe39ce4 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Mon, 9 Mar 2026 11:18:57 +0100 Subject: [PATCH 06/16] Ignore sonarqube error about cognitive complexity Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index d96e11f2..20d3f75a 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -260,7 +260,7 @@ def _check_quasi_newton_arguments(tol: float, maxiter: int) -> None: raise ValueError("maxiter must be greater than 0") -def quasi_newton( +def quasi_newton( # NOSONAR(S3776) func: Callable[[floatArrayLike], floatArrayLike], x0: floatArrayLike, tol: float = 1.48e-8, From 3d04327e8d5faf01b7af7be1270b42b62704bd7c Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:34:11 +0100 Subject: [PATCH 07/16] Refacto ConvectiveCooling and fix handling of wind_attack_angle - Make cigre ConvectiveCooling inherit ConvectiveCoolingBase like the others, so as to factor wind attack angle check and computations. - Handle missing case in wind attack angle computation (when it's nan or an array containing nans) - Enable orriding wind attack angle when initializing a solver - Rename _nu_forced -> _nusselt_forced and _nu_natural -> _nusselt_natural Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- .../power/cigre/convective_cooling.py | 31 +++++------- src/thermohl/power/convective_cooling.py | 47 +++++++++++++++---- src/thermohl/solver/parameters.py | 1 + .../test_power_cigre_convective_cooling.py | 32 ++++++------- .../docs/api-reference/parameters.md | 3 +- 5 files changed, 67 insertions(+), 47 deletions(-) diff --git a/src/thermohl/power/cigre/convective_cooling.py b/src/thermohl/power/cigre/convective_cooling.py index f2463634..a238519b 100644 --- a/src/thermohl/power/cigre/convective_cooling.py +++ b/src/thermohl/power/cigre/convective_cooling.py @@ -10,12 +10,13 @@ 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 -class ConvectiveCooling(PowerTerm): +class ConvectiveCooling(ConvectiveCoolingBase): """Convective cooling term.""" def __init__( @@ -46,6 +47,9 @@ def __init__( 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 @@ -53,20 +57,7 @@ def __init__( 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: - print( - "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: """ @@ -129,7 +120,7 @@ def _nu_forced( B1 * reynolds**n ) - def _nu_natural( + def _nusselt_natural( self, film_temperature: floatArrayLike, temperature_delta: floatArrayLike, @@ -188,8 +179,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 ( diff --git a/src/thermohl/power/convective_cooling.py b/src/thermohl/power/convective_cooling.py index 1451bb60..095dc87c 100644 --- a/src/thermohl/power/convective_cooling.py +++ b/src/thermohl/power/convective_cooling.py @@ -46,28 +46,55 @@ def __init__( wind_attack_angle: floatArrayLike = 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 + + def _check_arguments( + self, wind_azimuth: floatArrayLike, wind_attack_angle: floatArrayLike + ) -> None: + if ( + wind_attack_angle is None or np.isnan(wind_attack_angle).any() + ) 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: + if ( + wind_attack_angle is not None + and not np.isnan(wind_attack_angle).all() + and wind_azimuth is not None + ): print( "Warning: 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: + 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, diff --git a/src/thermohl/solver/parameters.py b/src/thermohl/solver/parameters.py index aeb91d83..2c7b27cb 100644 --- a/src/thermohl/solver/parameters.py +++ b/src/thermohl/solver/parameters.py @@ -50,6 +50,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 = 0 # nebulosity (1) self.albedo = 0.15 # albedo (1) # coefficient for air pollution from 0 (clean) to 1 (polluted) diff --git a/test/unit/power/cigre/test_power_cigre_convective_cooling.py b/test/unit/power/cigre/test_power_cigre_convective_cooling.py index 3a746f55..c3401112 100644 --- a/test/unit/power/cigre/test_power_cigre_convective_cooling.py +++ b/test/unit/power/cigre/test_power_cigre_convective_cooling.py @@ -45,12 +45,12 @@ conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_forced_float_value(convective_cooling, expected_type): +def test_nusselt_forced_float_value(convective_cooling, expected_type): Tf = 30.0 nu = 1.5e-5 expected_result = 17.3425 - result = convective_cooling._nu_forced(Tf, nu) + result = convective_cooling._nusselt_forced(Tf, nu) np.testing.assert_allclose(result, expected_result, rtol=1e-5) @@ -60,12 +60,12 @@ def test_nu_forced_float_value(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_forced_array_single_value(convective_cooling, expected_type): +def test_nusselt_forced_array_single_value(convective_cooling, expected_type): Tf = np.array([30.0]) nu = np.array([1.5e-5]) expected_result = np.array([17.3425]) - result = convective_cooling._nu_forced(Tf, nu) + result = convective_cooling._nusselt_forced(Tf, nu) np.testing.assert_allclose(result, expected_result, rtol=1e-5) @@ -75,12 +75,12 @@ def test_nu_forced_array_single_value(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_forced_array_values(convective_cooling, expected_type): +def test_nusselt_forced_array_values(convective_cooling, expected_type): Tf = np.array([30.0, 35.0]) nu = np.array([1.5e-5, 1.6e-5]) expected_result = np.array([17.3425, 16.6482]) - result = convective_cooling._nu_forced(Tf, nu) + result = convective_cooling._nusselt_forced(Tf, nu) np.testing.assert_allclose(result, expected_result, rtol=1e-5) @@ -90,14 +90,14 @@ def test_nu_forced_array_values(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_forced_boundary_conditions(convective_cooling, expected_type): +def test_nusselt_forced_boundary_conditions(convective_cooling, expected_type): Tf = np.array([30.0]) nu = np.array([1.5e-5]) convective_cooling.roughness_ratio = np.array([0.05]) convective_cooling.wind_speed = np.array([100.0]) expected_result = np.array([115.5214]) - result = convective_cooling._nu_forced(Tf, nu) + result = convective_cooling._nusselt_forced(Tf, nu) np.testing.assert_allclose(result, expected_result, rtol=1e-5) @@ -107,12 +107,12 @@ def test_nu_forced_boundary_conditions(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_natural_single_value(convective_cooling, expected_type): +def test_nusselt_natural_single_value(convective_cooling, expected_type): air_temperature = 50.0 Td = 25.0 nu = Air.kinematic_viscosity(air_temperature) - result = convective_cooling._nu_natural(air_temperature, Td, nu) + result = convective_cooling._nusselt_natural(air_temperature, Td, nu) assert isinstance(result, expected_type) np.testing.assert_allclose(result, [3.42404], rtol=1e-4) @@ -123,12 +123,12 @@ def test_nu_natural_single_value(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_natural_array_values(convective_cooling, expected_type): +def test_nusselt_natural_array_values(convective_cooling, expected_type): air_temperature = np.array([50.0, 60.0, 70.0]) Td = np.array([25.0, 35.0, 45.0]) nu = Air.kinematic_viscosity(air_temperature) - result = convective_cooling._nu_natural(air_temperature, Td, nu) + result = convective_cooling._nusselt_natural(air_temperature, Td, nu) assert isinstance(result, np.ndarray) assert result.shape == air_temperature.shape @@ -140,12 +140,12 @@ def test_nu_natural_array_values(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_natural_edge_case(convective_cooling, expected_type): +def test_nusselt_natural_edge_case(convective_cooling, expected_type): air_temperature = 0.0 Td = 0.0 nu = Air.kinematic_viscosity(air_temperature) - result = convective_cooling._nu_natural(air_temperature, Td, nu) + result = convective_cooling._nusselt_natural(air_temperature, Td, nu) assert isinstance(result, expected_type) assert result == 0.0 @@ -156,12 +156,12 @@ def test_nu_natural_edge_case(convective_cooling, expected_type): conv_cool_instances, ids=["ConvectiveCooling with arrays", "ConvectiveCooling with scalars"], ) -def test_nu_natural_high_values(convective_cooling, expected_type): +def test_nusselt_natural_high_values(convective_cooling, expected_type): air_temperature = 1000.0 Td = 500.0 nu = Air.kinematic_viscosity(air_temperature) - result = convective_cooling._nu_natural(air_temperature, Td, nu) + result = convective_cooling._nusselt_natural(air_temperature, Td, nu) assert isinstance(result, expected_type) np.testing.assert_allclose(result, [2.18853], rtol=1e-4) diff --git a/thermohl-docs/docs/api-reference/parameters.md b/thermohl-docs/docs/api-reference/parameters.md index 303740dc..f39d4320 100644 --- a/thermohl-docs/docs/api-reference/parameters.md +++ b/thermohl-docs/docs/api-reference/parameters.md @@ -32,7 +32,8 @@ units, default values and in which set of power terms they are used. | hour | 12 | N/A | yes | yes | yes | yes | hour of the day (float in[0, 24[) | | ambient_temperature | 15 | celsius | yes | yes | yes | yes | ambient temperature | | wind_speed | 0 | linear_mass.s⁻¹ | yes | yes | yes | yes | wind speed | -| wind_azimuth | 90 | degree | yes | yes | yes | yes | wind_azimuth (regarding north) | +| wind_azimuth | 90 | degree | yes | yes | yes | yes | wind azimuth (regarding north) | +| wind_attack_angle | None | radian | yes | yes | yes | yes | (absolute) angle between the wind and the cable. Must be comprised between 0 and π / 2. If provided, this overrides the angle computed from wind_azimuth and cable_azimuth. | | albedo | 0.8 | N/A | yes | no | no | no | albedo | | turbidity | 0.1 | N/A | no | yes | no | no | coefficient for air pollution from 0 (clean) to 1 (polluted) | | transit | 100 | A | yes | yes | yes | yes | transit intensity | From 8fd9602f717868150170c3baac58a31042892d88 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:34:22 +0100 Subject: [PATCH 08/16] Various fixes Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 28 +++++++++++++++++++++------- src/thermohl/utils.py | 2 +- test/unit/solver/test_slv1t.py | 18 +++++++++++++++--- thermohl-docs/docs/user-guide.md | 2 +- 4 files changed, 38 insertions(+), 12 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index eca9d131..685487ac 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -279,20 +279,23 @@ def reduced_intensity( and a faulty sleeve. Args: - measured_temperature_difference (float): The measured temperature difference between the cable surface and the sleeve. - measured_intensity (float): The measuredintensity at which the temperature difference was measured. - ambient_temperature (Optional[float]): The ambient temperature. Default is 30. - wind_speed (Optional[float]): The wind speed. Default is 0.6. - measured_solar_irradiance (Optional[float]): The measured solar irradiance. Default is 600. - max_conductor_temperature (Optional[float]): The maximum conductor temperature. Default is 100. + 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. + measured_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 args that will be modified so as to be able to restore them at the end of the computation solver_transit = self.args.transit solver_ambient_temperature = self.args.ambient_temperature solver_wind_speed = self.args.wind_speed solver_measured_solar_irradiance = self.args.measured_solar_irradiance + solver_has_wind_attack_angle = hasattr(self.args, "wind_attack_angle") + if solver_has_wind_attack_angle: + solver_wind_attack_angle = self.args.wind_attack_angle - # Set default values for reduced intensity computation. + # Set args default values for reduced intensity computation. # These differ from those used for the other computations. self._set_default_reduced_intensity_args( ambient_temperature, wind_speed, measured_solar_irradiance @@ -302,6 +305,9 @@ def reduced_intensity( if max_conductor_temperature is None: max_conductor_temperature = np.ones_like(measured_intensity) * 100.0 + self.args.wind_attack_angle = 90.0 + self.convective_cooling.__init__(**self.args.__dict__) + def conductor_temperature(transit): self.args.transit = transit self.joule_heating.__init__(**self.args.__dict__) @@ -324,8 +330,16 @@ def f(transit): # Restore previous args self.args.transit = solver_transit + # Update joule heating with restored transit + self.joule_heating.__init__(**self.args.__dict__) self.args.ambient_temperature = solver_ambient_temperature self.args.wind_speed = solver_wind_speed self.args.measured_solar_irradiance = solver_measured_solar_irradiance + if solver_has_wind_attack_angle: + self.args.wind_attack_angle = solver_wind_attack_angle + elif hasattr(self.args, "wind_attack_angle"): + del self.args.wind_attack_angle + # Update convective cooling with restored wind_attack_angle + self.convective_cooling.__init__(**self.args.__dict__) return reduced_intensity diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index 20d3f75a..541ab8f6 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -197,7 +197,7 @@ def _array_quasi_newton( A vectorized version of secant method for arrays. Do not use this method directly. This method is called from `quasi_newton` - when ``np.size(x0) > 1`` is ``True``. For docstring, see `quasi_newton`. + when ``np.size(x0) > 1`` is ``True``. Heavily inspired by the implementation of the SciPy library. """ diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 540e31cd..f9922405 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -220,10 +220,10 @@ def test_steady_intensity_custom_params(solver): def test_reduced_intensity_scalar(): args = { "max_len": lambda: 1, - VariableType.TRANSIT.value: 0, + VariableType.TRANSIT.value: 20, "ambient_temperature": 25, "wind_speed": 0, - "wind_azimuth": 0, + "wind_azimuth": 30.0, "ambient_pressure": 101325, "relative_humidity": 50, "precipitation_rate": 0, @@ -253,6 +253,10 @@ def test_reduced_intensity_scalar(): args.update(cable) solver = create_solver(args) + + initial_wind_attack_angle = solver.convective_cooling.wind_attack_angle + initial_transit = 20 + result = solver.reduced_intensity( measured_temperature_difference=10.0, measured_intensity=360.0, @@ -260,11 +264,15 @@ def test_reduced_intensity_scalar(): assert isinstance(result, np.float64) - # Check that solver args have not been changed + # Check that solver args and power term attributes have not been changed assert np.isnan(solver.args.measured_solar_irradiance) assert solver.args.ambient_temperature == args["ambient_temperature"] assert solver.args.wind_speed == args["wind_speed"] assert solver.args.transit == args["transit"] + assert np.isclose(solver.args.wind_attack_angle, np.deg2rad(30.0)) + + assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle + assert solver.joule_heating.transit == initial_transit def test_reduced_intensity_array(): @@ -274,6 +282,7 @@ def test_reduced_intensity_array(): "ambient_temperature": np.array([25, 25]), "wind_speed": np.array([0, 0]), "wind_azimuth": np.array([0, 0]), + "wind_attack_angle": np.array([45.0, 45.0]), "ambient_pressure": np.array([101325, 101325]), "relative_humidity": np.array([50, 50]), "precipitation_rate": np.array([0, 0]), @@ -313,3 +322,6 @@ def test_reduced_intensity_array(): assert len(result) == 2 assert not np.isnan(result[0]) assert not np.isnan(result[1]) + + # Assert solver wind attack angle has not been changed + assert np.allclose(solver.args.wind_attack_angle, args["wind_attack_angle"]) diff --git a/thermohl-docs/docs/user-guide.md b/thermohl-docs/docs/user-guide.md index 90965000..a2cf49d9 100644 --- a/thermohl-docs/docs/user-guide.md +++ b/thermohl-docs/docs/user-guide.md @@ -161,7 +161,7 @@ However, if a sleeve is faulty, it can overheat. ThermoHL enables the user to compute a reduced intensity limit (called reduced IST) so that the interface between the sleeve and the cable doesn't exceed a given maximum temperature. Input data are -- the mesaured temperature difference between the hotspot (at the interface between sleeve and cable) +- the measured temperature difference between the hotspot (at the interface between sleeve and cable) and the cable far away from the sleeve - the measured transit. From c9ce67ed9e44fd951f3d3a4c93446903be319e3f Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Mon, 9 Mar 2026 16:47:50 +0100 Subject: [PATCH 09/16] Remove dead code and add test Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 9 ++--- test/unit/solver/test_slv1t.py | 64 ++++++++++++++++++++++++++++++++-- 2 files changed, 64 insertions(+), 9 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 685487ac..5caec032 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -291,9 +291,7 @@ def reduced_intensity( solver_ambient_temperature = self.args.ambient_temperature solver_wind_speed = self.args.wind_speed solver_measured_solar_irradiance = self.args.measured_solar_irradiance - solver_has_wind_attack_angle = hasattr(self.args, "wind_attack_angle") - if solver_has_wind_attack_angle: - solver_wind_attack_angle = self.args.wind_attack_angle + solver_wind_attack_angle = self.args.wind_attack_angle # Set args default values for reduced intensity computation. # These differ from those used for the other computations. @@ -335,10 +333,7 @@ def f(transit): self.args.ambient_temperature = solver_ambient_temperature self.args.wind_speed = solver_wind_speed self.args.measured_solar_irradiance = solver_measured_solar_irradiance - if solver_has_wind_attack_angle: - self.args.wind_attack_angle = solver_wind_attack_angle - elif hasattr(self.args, "wind_attack_angle"): - del self.args.wind_attack_angle + self.args.wind_attack_angle = solver_wind_attack_angle # Update convective cooling with restored wind_attack_angle self.convective_cooling.__init__(**self.args.__dict__) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index f9922405..103f7207 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -217,7 +217,7 @@ def test_steady_intensity_custom_params(solver): assert VariableType.TRANSIT in result.columns -def test_reduced_intensity_scalar(): +def test_reduced_intensity_scalar_using_default_args(): args = { "max_len": lambda: 1, VariableType.TRANSIT.value: 20, @@ -251,15 +251,75 @@ def test_reduced_intensity_scalar(): "emissivity": 0.8, } args.update(cable) + solver = create_solver(args) + + initial_wind_attack_angle = solver.convective_cooling.wind_attack_angle + initial_transit = args[VariableType.TRANSIT.value] + + result = solver.reduced_intensity( + measured_temperature_difference=10.0, + measured_intensity=360.0, + ) + + assert isinstance(result, np.float64) + # Check that solver args and power term attributes have not been changed + assert np.isnan(solver.args.measured_solar_irradiance) + assert solver.args.ambient_temperature == args["ambient_temperature"] + assert solver.args.wind_speed == args["wind_speed"] + assert solver.args.transit == args["transit"] + assert np.isclose(solver.args.wind_attack_angle, np.deg2rad(30.0)) + + assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle + assert solver.joule_heating.transit == initial_transit + + +def test_reduced_intensity_scalar_providing_custom_args(): + args = { + "max_len": lambda: 1, + VariableType.TRANSIT.value: 20, + "ambient_temperature": 25, + "wind_speed": 0, + "wind_azimuth": 30.0, + "ambient_pressure": 101325, + "relative_humidity": 50, + "precipitation_rate": 0, + "month": 1, + "day": 1, + "hour": 0, + "latitude": 48.0, + "longitude": 2.3, + "altitude": 50.0, + "cable_azimuth": 0.0, + "solar_absorptivity": 0.5, + } + cable = { + "linear_mass": 1.0, + "heat_capacity": 1.0, + "outer_diameter": 3.186e-02, + "core_diameter": 0.000, + "outer_area": 6.004e-4, + "core_area": 0.000, + "magnetic_coeff": 1.000, + "magnetic_coeff_per_a": 0.000, + "temperature_coeff_linear": 3.600e-3, + "temperature_coeff_quadratic": 8.000e-7, + "linear_resistance_dc_20c": 5.540e-5, + "emissivity": 0.8, + } + args.update(cable) solver = create_solver(args) initial_wind_attack_angle = solver.convective_cooling.wind_attack_angle - initial_transit = 20 + initial_transit = args[VariableType.TRANSIT.value] result = solver.reduced_intensity( measured_temperature_difference=10.0, measured_intensity=360.0, + ambient_temperature=25.0, + wind_speed=4.0, + measured_solar_irradiance=800.0, + max_conductor_temperature=120.0, ) assert isinstance(result, np.float64) From 46c4aa82b94926b2566d3625dd8149fb9612b4b2 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Wed, 11 Mar 2026 09:23:37 +0100 Subject: [PATCH 10/16] Add tests Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/utils.py | 6 ++- test/unit/test_utils.py | 94 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 96 insertions(+), 4 deletions(-) diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index 541ab8f6..a73c488b 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -253,11 +253,13 @@ def _array_quasi_newton( return p -def _check_quasi_newton_arguments(tol: float, maxiter: int) -> None: +def _check_quasi_newton_arguments(tol: float, maxiter: int, rtol: float) -> None: if tol <= 0: raise ValueError(f"tol too small ({tol:g} <= 0)") if maxiter < 1: raise ValueError("maxiter must be greater than 0") + if rtol < 0: + raise ValueError(f"rtol too small ({rtol:g} < 0)") def quasi_newton( # NOSONAR(S3776) @@ -271,7 +273,7 @@ def quasi_newton( # NOSONAR(S3776) Heavily inspired by the implementation of optimize.newton in the SciPy library. """ - _check_quasi_newton_arguments(tol, maxiter) + _check_quasi_newton_arguments(tol, maxiter, rtol) if np.size(x0) > 1: return _array_quasi_newton(func, x0, tol, maxiter) diff --git a/test/unit/test_utils.py b/test/unit/test_utils.py index ae648fe5..26ec3bb5 100644 --- a/test/unit/test_utils.py +++ b/test/unit/test_utils.py @@ -8,7 +8,7 @@ import numpy as np import pytest -from thermohl.utils import bisect_v, quasi_newton_2d +from thermohl.utils import bisect_v, quasi_newton, quasi_newton_2d _nprs = 3141592654 @@ -102,7 +102,97 @@ def f(x): ) -# +def test_quasi_newton_wrong_tol(): + def f(x): + return x**2 - 2 + + with pytest.raises(ValueError): + quasi_newton( + f, + x0=1.0, + tol=-1.0e-6, # Invalid tolerance + ) + + +def test_quasi_newton_wrong_maxiter(): + def f(x): + return x**2 - 2 + + with pytest.raises(ValueError): + quasi_newton( + f, + x0=1.0, + maxiter=0, # Invalid max iterations + ) + + +def test_quasi_newton_wrong_rtol(): + def f(x): + return x**2 - 2 + + with pytest.raises(ValueError): + quasi_newton( + f, + x0=1.0, + rtol=-1.0e-6, # Invalid relative tolerance + ) + + +def test_quasi_newton_scalar_increasing(): + def f(x): + return np.log(x) + + root = quasi_newton( + f, + x0=0.5, + ) + assert np.isclose(root, 1.0) + + +def test_quasi_newton_scalar_decreasing(): + def f(x): + return -np.log(x) + + root = quasi_newton( + f, + x0=0.5, + ) + assert np.isclose(root, 1.0) + + +def test_quasi_newton_scalar_no_convergence(): + def f(x): + return x**2 + 2 + + with pytest.raises(RuntimeError): + quasi_newton( + f, + x0=1.0, + ) + + +def test_quasi_newton_array_no_convergence(): + def f(x: np.ndarray) -> np.ndarray: + return x**2 + np.array([1, 2]) + + with pytest.raises(RuntimeError): + quasi_newton( + f, + x0=np.array([1.0, 1.0]), + ) + + +def test_quasi_newton_array_mixed(): + def f(x: np.ndarray) -> np.ndarray: + return x**2 + np.array([-1, 2]) + + result = quasi_newton( + f, + x0=np.array([1.0, 1.0]), + ) + assert np.isclose(result[0], 1.0) + + def test_quasi_newton_2d_convergence(): np.random.seed(_nprs) size = 10 From ef183763e9ff023a54f2fd7880639c0bc115579b Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Thu, 19 Mar 2026 08:51:12 +0100 Subject: [PATCH 11/16] Add type hints Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 2 +- test/unit/solver/test_slv1t.py | 6 +++--- test/unit/test_utils.py | 16 ++++++++-------- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 5caec032..efc6b871 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -272,7 +272,7 @@ def reduced_intensity( wind_speed: Optional[floatArrayLike] = None, measured_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 diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 103f7207..747ea239 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -217,7 +217,7 @@ def test_steady_intensity_custom_params(solver): assert VariableType.TRANSIT in result.columns -def test_reduced_intensity_scalar_using_default_args(): +def test_reduced_intensity_scalar_using_default_args() -> None: args = { "max_len": lambda: 1, VariableType.TRANSIT.value: 20, @@ -274,7 +274,7 @@ def test_reduced_intensity_scalar_using_default_args(): assert solver.joule_heating.transit == initial_transit -def test_reduced_intensity_scalar_providing_custom_args(): +def test_reduced_intensity_scalar_providing_custom_args() -> None: args = { "max_len": lambda: 1, VariableType.TRANSIT.value: 20, @@ -335,7 +335,7 @@ def test_reduced_intensity_scalar_providing_custom_args(): assert solver.joule_heating.transit == initial_transit -def test_reduced_intensity_array(): +def test_reduced_intensity_array() -> None: args = { "max_len": lambda: 2, VariableType.TRANSIT.value: np.array([0, 0]), diff --git a/test/unit/test_utils.py b/test/unit/test_utils.py index 26ec3bb5..57df1da0 100644 --- a/test/unit/test_utils.py +++ b/test/unit/test_utils.py @@ -102,7 +102,7 @@ def f(x): ) -def test_quasi_newton_wrong_tol(): +def test_quasi_newton_wrong_tol() -> None: def f(x): return x**2 - 2 @@ -114,7 +114,7 @@ def f(x): ) -def test_quasi_newton_wrong_maxiter(): +def test_quasi_newton_wrong_maxiter() -> None: def f(x): return x**2 - 2 @@ -126,7 +126,7 @@ def f(x): ) -def test_quasi_newton_wrong_rtol(): +def test_quasi_newton_wrong_rtol() -> None: def f(x): return x**2 - 2 @@ -138,7 +138,7 @@ def f(x): ) -def test_quasi_newton_scalar_increasing(): +def test_quasi_newton_scalar_increasing() -> None: def f(x): return np.log(x) @@ -149,7 +149,7 @@ def f(x): assert np.isclose(root, 1.0) -def test_quasi_newton_scalar_decreasing(): +def test_quasi_newton_scalar_decreasing() -> None: def f(x): return -np.log(x) @@ -160,7 +160,7 @@ def f(x): assert np.isclose(root, 1.0) -def test_quasi_newton_scalar_no_convergence(): +def test_quasi_newton_scalar_no_convergence() -> None: def f(x): return x**2 + 2 @@ -171,7 +171,7 @@ def f(x): ) -def test_quasi_newton_array_no_convergence(): +def test_quasi_newton_array_no_convergence() -> None: def f(x: np.ndarray) -> np.ndarray: return x**2 + np.array([1, 2]) @@ -182,7 +182,7 @@ def f(x: np.ndarray) -> np.ndarray: ) -def test_quasi_newton_array_mixed(): +def test_quasi_newton_array_mixed() -> None: def f(x: np.ndarray) -> np.ndarray: return x**2 + np.array([-1, 2]) From 12723112ddd189b67a80a1ea756619e6c58400be Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Fri, 10 Apr 2026 10:48:04 +0200 Subject: [PATCH 12/16] Formatting + improvements after review Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 22 +++++++++---------- test/unit/solver/test_slv1t.py | 6 ++--- .../docs/examples/ex_ist_reduction.ipynb | 12 +++++----- thermohl-docs/docs/user-guide.md | 2 +- 4 files changed, 22 insertions(+), 20 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index efc6b871..1ab52701 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -16,7 +16,7 @@ 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 bisect_v, quasi_newton +from thermohl.utils import quasi_newton class Solver1T(Solver_): @@ -243,7 +243,7 @@ def _set_default_reduced_intensity_args( self, ambient_temperature: Optional[floatArrayLike], wind_speed: Optional[floatArrayLike], - measured_solar_irradiance: Optional[floatArrayLike], + measured_global_radiation: Optional[floatArrayLike], ): if ambient_temperature is None: print( @@ -257,12 +257,12 @@ def _set_default_reduced_intensity_args( wind_speed = 0.6 self.args.wind_speed = wind_speed - if measured_solar_irradiance is None: + if measured_global_radiation is None: print( - "WARNING: measured_solar_irradiance is not set. Using default value of 600 W/m²." + "WARNING: measured_global_radiation is not set. Using default value of 600 W/m²." ) - measured_solar_irradiance = 600.0 - self.args.measured_solar_irradiance = measured_solar_irradiance + measured_global_radiation = 600.0 + self.args.measured_global_radiation = measured_global_radiation def reduced_intensity( self, @@ -270,7 +270,7 @@ def reduced_intensity( measured_intensity: floatArrayLike, ambient_temperature: Optional[floatArrayLike] = None, wind_speed: Optional[floatArrayLike] = None, - measured_solar_irradiance: Optional[floatArrayLike] = None, + measured_global_radiation: Optional[floatArrayLike] = None, max_conductor_temperature: Optional[floatArrayLike] = None, ) -> floatArrayLike: """ @@ -283,20 +283,20 @@ def reduced_intensity( 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. - measured_solar_irradiance (Optional[float | np.ndarray]): The measured solar irradiance. Default is 600. + measured_global_radiation (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 args that will be modified so as to be able to restore them at the end of the computation solver_transit = self.args.transit solver_ambient_temperature = self.args.ambient_temperature solver_wind_speed = self.args.wind_speed - solver_measured_solar_irradiance = self.args.measured_solar_irradiance + solver_measured_solar_irradiance = self.args.measured_global_radiation solver_wind_attack_angle = self.args.wind_attack_angle # Set args default values for reduced intensity computation. # These differ from those used for the other computations. self._set_default_reduced_intensity_args( - ambient_temperature, wind_speed, measured_solar_irradiance + ambient_temperature, wind_speed, measured_global_radiation ) # Set default value for max_conductor_temperature if not provided. @@ -332,7 +332,7 @@ def f(transit): self.joule_heating.__init__(**self.args.__dict__) self.args.ambient_temperature = solver_ambient_temperature self.args.wind_speed = solver_wind_speed - self.args.measured_solar_irradiance = solver_measured_solar_irradiance + self.args.measured_global_radiation = solver_measured_solar_irradiance self.args.wind_attack_angle = solver_wind_attack_angle # Update convective cooling with restored wind_attack_angle self.convective_cooling.__init__(**self.args.__dict__) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 747ea239..ad1f3351 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -264,7 +264,7 @@ def test_reduced_intensity_scalar_using_default_args() -> None: assert isinstance(result, np.float64) # Check that solver args and power term attributes have not been changed - assert np.isnan(solver.args.measured_solar_irradiance) + assert np.isnan(solver.args.measured_global_radiation) assert solver.args.ambient_temperature == args["ambient_temperature"] assert solver.args.wind_speed == args["wind_speed"] assert solver.args.transit == args["transit"] @@ -318,14 +318,14 @@ def test_reduced_intensity_scalar_providing_custom_args() -> None: measured_intensity=360.0, ambient_temperature=25.0, wind_speed=4.0, - measured_solar_irradiance=800.0, + measured_global_radiation=800.0, max_conductor_temperature=120.0, ) assert isinstance(result, np.float64) # Check that solver args and power term attributes have not been changed - assert np.isnan(solver.args.measured_solar_irradiance) + assert np.isnan(solver.args.measured_global_radiation) assert solver.args.ambient_temperature == args["ambient_temperature"] assert solver.args.wind_speed == args["wind_speed"] assert solver.args.transit == args["transit"] diff --git a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb index 789a16af..af35c738 100644 --- a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb +++ b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb @@ -10,8 +10,8 @@ "import numpy as np\n", "\n", "from thermohl import solver\n", - "from thermohl.solver.enums.heat_equation_type import HeatEquationType\n", - "from thermohl.solver.enums.variable_type import VariableType\n" + "from thermohl.solver import HeatEquationType\n", + "from thermohl.solver.entities import VariableType\n" ] }, { @@ -50,13 +50,15 @@ "text": [ "WARNING: ambient_temperature is not set. Using default value of 30 °C.\n", "WARNING: wind_speed is not set. Using default value of 0.6 m/s.\n", - "WARNING: measured_solar_irradiance is not set. Using default value of 600 W/m².\n" + "WARNING: measured_global_radiation is not set. Using default value of 600 W/m².\n", + "Warning: both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored.\n", + "Warning: both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored.\n" ] }, { "data": { "text/plain": [ - "np.float64(830.1154202065042)" + "np.float64(850.0270802708803)" ] }, "execution_count": 3, @@ -67,7 +69,7 @@ "source": [ "slvr = solver.ieee(\n", " args,\n", - " heat_equation=HeatEquationType.WITH_ONE_TEMPERATURE,\n", + " heat_equation=HeatEquationType.ONE_TEMPERATURE,\n", ")\n", "\n", "slvr.reduced_intensity(\n", diff --git a/thermohl-docs/docs/user-guide.md b/thermohl-docs/docs/user-guide.md index a2cf49d9..65aaad8c 100644 --- a/thermohl-docs/docs/user-guide.md +++ b/thermohl-docs/docs/user-guide.md @@ -174,6 +174,6 @@ computations: |----|----| | ambient_temperature | 30.0 °C | | wind_speed | 0.6 m/s | -| measured_solar_irradiance | 600.0 W/m² | +| measured_global_radiation | 600.0 W/m² | TODO: precision on wind_speed and wind_angle if needed From 74a55d80cf1fd8a7afb734069220d998593421e9 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Fri, 10 Apr 2026 15:18:07 +0200 Subject: [PATCH 13/16] Use logger instead of print Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/solver/slv1t.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 1ab52701..e249809d 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -5,6 +5,7 @@ # 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 Dict, Any, Optional @@ -19,6 +20,9 @@ from thermohl.utils import quasi_newton +logger = logging.getLogger(__name__) + + class Solver1T(Solver_): def steady_temperature( self, @@ -246,20 +250,20 @@ def _set_default_reduced_intensity_args( measured_global_radiation: Optional[floatArrayLike], ): if ambient_temperature is None: - print( - "WARNING: ambient_temperature is not set. Using default value of 30 °C." + 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: - print("WARNING: wind_speed is not set. Using default value of 0.6 m/s.") + 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 measured_global_radiation is None: - print( - "WARNING: measured_global_radiation is not set. Using default value of 600 W/m²." + logger.warning( + "measured_global_radiation is not set. Using default value of 600 W/m²." ) measured_global_radiation = 600.0 self.args.measured_global_radiation = measured_global_radiation From 0f5a911ec3d97112429fe4a5c95bf4a07c59db2f Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Tue, 14 Apr 2026 09:57:05 +0200 Subject: [PATCH 14/16] Improvements after review Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- .../power/cigre/convective_cooling.py | 7 +++-- src/thermohl/power/convective_cooling.py | 30 +++++++++++++++---- src/thermohl/solver/slv1t.py | 4 +-- 3 files changed, 30 insertions(+), 11 deletions(-) diff --git a/src/thermohl/power/cigre/convective_cooling.py b/src/thermohl/power/cigre/convective_cooling.py index e5ca22b0..563db3c0 100644 --- a/src/thermohl/power/cigre/convective_cooling.py +++ b/src/thermohl/power/cigre/convective_cooling.py @@ -31,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, ): @@ -45,7 +45,8 @@ 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. diff --git a/src/thermohl/power/convective_cooling.py b/src/thermohl/power/convective_cooling.py index 8d5bdb08..784296f9 100644 --- a/src/thermohl/power/convective_cooling.py +++ b/src/thermohl/power/convective_cooling.py @@ -46,8 +46,8 @@ 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) @@ -63,13 +63,17 @@ def __init__( self.dynamic_viscosity = dynamic_viscosity self.thermal_conductivity = thermal_conductivity + @classmethod def _check_arguments( - self, wind_azimuth: floatArrayLike, wind_attack_angle: floatArrayLike + 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: + 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 not np.isnan(wind_attack_angle).all() @@ -87,6 +91,20 @@ def _set_wind_attack_angle( ) -> 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( diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index e249809d..b2c82534 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -305,7 +305,7 @@ def reduced_intensity( # Set default value for max_conductor_temperature if not provided. if max_conductor_temperature is None: - max_conductor_temperature = np.ones_like(measured_intensity) * 100.0 + max_conductor_temperature = np.full_like(measured_intensity, 100.0) self.args.wind_attack_angle = 90.0 self.convective_cooling.__init__(**self.args.__dict__) @@ -326,7 +326,7 @@ def sleeve_temperature(transit): def f(transit): return sleeve_temperature(transit) - max_conductor_temperature - x0 = np.ones_like(measured_intensity) * 100 + x0 = np.full_like(measured_intensity, 100.0) reduced_intensity = quasi_newton(f, x0=x0) From 8c78c26433f56415d68178e561c5b88de13c9dc6 Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Thu, 16 Apr 2026 08:57:57 +0200 Subject: [PATCH 15/16] notebooks: Update Ist reduction notebook Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- .../docs/examples/ex_ist_reduction.ipynb | 36 ++++--------------- 1 file changed, 7 insertions(+), 29 deletions(-) diff --git a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb index af35c738..dc8ffb5c 100644 --- a/thermohl-docs/docs/examples/ex_ist_reduction.ipynb +++ b/thermohl-docs/docs/examples/ex_ist_reduction.ipynb @@ -2,11 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "8c937b14", "metadata": {}, "outputs": [], "source": [ + "from datetime import datetime\n", + "\n", "import numpy as np\n", "\n", "from thermohl import solver\n", @@ -16,7 +18,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "81557e4b", "metadata": {}, "outputs": [], @@ -32,40 +34,16 @@ " \"precipitation_rate\": np.array([0]),\n", " \"linear_mass\": 1.0,\n", " \"heat_capacity\": 1.0,\n", - " \"month\": 1,\n", - " \"day\": 1,\n", - " \"hour\": 0,\n", + " \"datetime_utc\": datetime(2026, 1, 1, 0, 0, 0),\n", "}" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "0cf8da43", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "WARNING: ambient_temperature is not set. Using default value of 30 °C.\n", - "WARNING: wind_speed is not set. Using default value of 0.6 m/s.\n", - "WARNING: measured_global_radiation is not set. Using default value of 600 W/m².\n", - "Warning: both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored.\n", - "Warning: both wind_attack_angle and wind_azimuth are provided. wind_azimuth will be ignored.\n" - ] - }, - { - "data": { - "text/plain": [ - "np.float64(850.0270802708803)" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "slvr = solver.ieee(\n", " args,\n", From f1de9d61dcfa221606d835b9b4704485497cd75a Mon Sep 17 00:00:00 2001 From: ai-qui <184963772+ai-qui@users.noreply.github.com> Date: Wed, 22 Apr 2026 14:03:47 +0200 Subject: [PATCH 16/16] Fix handling of solar heating/solar irradiance The user wants the be able to set the solar irradiance, not the measured global radiation. Adds new solar heating power term `FixedSolarIrradianceSolarHeating` to be able to set the solar irradiance instead of computing it. Signed-off-by: ai-qui <184963772+ai-qui@users.noreply.github.com> --- src/thermohl/power/__init__.py | 3 +- src/thermohl/power/solar_heating.py | 18 ++++++++++ src/thermohl/solver/slv1t.py | 51 +++++++++++++++++------------ test/unit/solver/test_slv1t.py | 5 ++- 4 files changed, 54 insertions(+), 23 deletions(-) diff --git a/src/thermohl/power/__init__.py b/src/thermohl/power/__init__.py index 3428bf56..8391ba2c 100644 --- a/src/thermohl/power/__init__.py +++ b/src/thermohl/power/__init__.py @@ -8,7 +8,7 @@ from .radiative_cooling import RadiativeCoolingBase from .power_term import PowerTerm -from .solar_heating import _SRad, SolarHeatingBase +from .solar_heating import _SRad, SolarHeatingBase, FixedSolarIrradianceSolarHeating __all__ = [ @@ -16,4 +16,5 @@ "PowerTerm", "_SRad", "SolarHeatingBase", + "FixedSolarIrradianceSolarHeating", ] diff --git a/src/thermohl/power/solar_heating.py b/src/thermohl/power/solar_heating.py index e7d0bfab..06d00a5a 100644 --- a/src/thermohl/power/solar_heating.py +++ b/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 diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index b2c82534..64116868 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -13,6 +13,7 @@ import pandas as pd 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 @@ -247,7 +248,7 @@ def _set_default_reduced_intensity_args( self, ambient_temperature: Optional[floatArrayLike], wind_speed: Optional[floatArrayLike], - measured_global_radiation: Optional[floatArrayLike], + solar_irradiance: Optional[floatArrayLike], ): if ambient_temperature is None: logger.warning( @@ -261,12 +262,12 @@ def _set_default_reduced_intensity_args( wind_speed = 0.6 self.args.wind_speed = wind_speed - if measured_global_radiation is None: + if solar_irradiance is None: logger.warning( - "measured_global_radiation is not set. Using default value of 600 W/m²." + "solar_irradiance is not set. Using default value of 600 W/m²." ) - measured_global_radiation = 600.0 - self.args.measured_global_radiation = measured_global_radiation + solar_irradiance = 600.0 + return solar_irradiance def reduced_intensity( self, @@ -274,7 +275,7 @@ def reduced_intensity( measured_intensity: floatArrayLike, ambient_temperature: Optional[floatArrayLike] = None, wind_speed: Optional[floatArrayLike] = None, - measured_global_radiation: Optional[floatArrayLike] = None, + solar_irradiance: Optional[floatArrayLike] = None, max_conductor_temperature: Optional[floatArrayLike] = None, ) -> floatArrayLike: """ @@ -287,20 +288,22 @@ def reduced_intensity( 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. - measured_global_radiation (Optional[float | np.ndarray]): The measured solar irradiance. Default is 600. + 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 args that will be modified so as to be able to restore them at the end of the computation - solver_transit = self.args.transit - solver_ambient_temperature = self.args.ambient_temperature - solver_wind_speed = self.args.wind_speed - solver_measured_solar_irradiance = self.args.measured_global_radiation - solver_wind_attack_angle = self.args.wind_attack_angle + # 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. - self._set_default_reduced_intensity_args( - ambient_temperature, wind_speed, measured_global_radiation + solar_irradiance = self._set_default_reduced_intensity_args( + ambient_temperature, wind_speed, solar_irradiance ) # Set default value for max_conductor_temperature if not provided. @@ -310,6 +313,12 @@ def reduced_intensity( 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__) @@ -330,15 +339,15 @@ def f(transit): reduced_intensity = quasi_newton(f, x0=x0) - # Restore previous args - self.args.transit = solver_transit + # 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 = solver_ambient_temperature - self.args.wind_speed = solver_wind_speed - self.args.measured_global_radiation = solver_measured_solar_irradiance - self.args.wind_attack_angle = solver_wind_attack_angle + 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 diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index ad1f3351..96617088 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -12,6 +12,7 @@ from thermohl import power from thermohl.solver.entities import PowerType, VariableType +from thermohl.power.rte.solar_heating import SolarHeating as RteSolarHeating from thermohl.solver.slv1t import Solver1T @@ -272,6 +273,7 @@ def test_reduced_intensity_scalar_using_default_args() -> None: assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle assert solver.joule_heating.transit == initial_transit + assert isinstance(solver.solar_heating, RteSolarHeating) def test_reduced_intensity_scalar_providing_custom_args() -> None: @@ -318,7 +320,7 @@ def test_reduced_intensity_scalar_providing_custom_args() -> None: measured_intensity=360.0, ambient_temperature=25.0, wind_speed=4.0, - measured_global_radiation=800.0, + solar_irradiance=800.0, max_conductor_temperature=120.0, ) @@ -333,6 +335,7 @@ def test_reduced_intensity_scalar_providing_custom_args() -> None: assert solver.convective_cooling.wind_attack_angle == initial_wind_attack_angle assert solver.joule_heating.transit == initial_transit + assert isinstance(solver.solar_heating, RteSolarHeating) def test_reduced_intensity_array() -> None: