diff --git a/linerate/equations/heat_capacity.py b/linerate/equations/heat_capacity.py new file mode 100644 index 0000000..f93581f --- /dev/null +++ b/linerate/equations/heat_capacity.py @@ -0,0 +1,20 @@ +from linerate.units import ( + Celsius, + JoulePerKilogramPerKelvin, + KilogramPerMeter, + PerKelvin, + JoulePerMeter, +) + + +def calculate_heat_capacity_per_unit_length( + mass_per_unit_length: KilogramPerMeter, + specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin, + specific_heat_capacity_coefficient: PerKelvin, + conductor_temperature: Celsius, +) -> JoulePerMeter: + m = mass_per_unit_length + c = specific_heat_capacity_at_20_celsius + T = conductor_temperature + beta = specific_heat_capacity_coefficient + return m * c * (1 + (T - 20) * beta) diff --git a/linerate/models/cigre601.py b/linerate/models/cigre601.py index 7ab9c6e..7fbcc1b 100644 --- a/linerate/models/cigre601.py +++ b/linerate/models/cigre601.py @@ -136,9 +136,8 @@ def compute_temperature_gradient( n = self.span.num_conductors T_c = conductor_temperature I = current / n # noqa - R = self.compute_resistance(conductor_temperature=T_c, current=I) return cigre601.convective_cooling.compute_temperature_gradient( - total_heat_gain=I * R, + total_heat_gain=self.compute_joule_heating(conductor_temperature=T_c, current=I), conductor_thermal_conductivity=self.span.conductor.thermal_conductivity, # type: ignore # noqa core_diameter=self.span.conductor.core_diameter, conductor_diameter=self.span.conductor.conductor_diameter, diff --git a/linerate/models/thermal_model.py b/linerate/models/thermal_model.py index 61bab94..62297b0 100644 --- a/linerate/models/thermal_model.py +++ b/linerate/models/thermal_model.py @@ -1,10 +1,13 @@ from abc import ABC, abstractmethod from typing import Dict +import numpy as np +from functools import partial from linerate import solver -from linerate.equations import joule_heating, radiative_cooling -from linerate.types import BaseWeather, Span -from linerate.units import Ampere, Celsius, OhmPerMeter, WattPerMeter +from linerate.equations import joule_heating, radiative_cooling, heat_capacity +from linerate.solver import solve_ivp_forward_euler +from linerate.types import BaseWeather, ConductorWithTransientData, Span +from linerate.units import Ampere, Celsius, Duration, OhmPerMeter, WattPerMeter, JoulePerKilogramPerKelvin def _copy_method_docstring(parent_class): @@ -269,3 +272,126 @@ def compute_conductor_temperature( tolerance=tolerance, ) return T + + def compute_heat_capacity_per_unit_length(self, conductor_temperature: Celsius) -> JoulePerKilogramPerKelvin: + r""" + Compute the heat capacity of the conductor per unit length. + Parameters + ---------- + conductor_temperature: + :math:`T~\left[^\circ\text{C}\right]` Conductor temperature. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`c~\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\text{m}^{-1}\right]`. Conductor heat capacity per unit length. + + """ + conductor = self.span.conductor + if not isinstance(conductor, ConductorWithTransientData): + raise ValueError("Heat capacity data must be defined to calculate heat capacity.") + mc_s = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.steel_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.steel_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.steel_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + mc_a = heat_capacity.calculate_heat_capacity_per_unit_length( + mass_per_unit_length=conductor.aluminum_mass_per_unit_length, + specific_heat_capacity_at_20_celsius=conductor.aluminum_specific_heat_capacity_at_20_celsius, + specific_heat_capacity_coefficient=conductor.aluminum_specific_heat_capacity_temperature_coefficient, + conductor_temperature=conductor_temperature, + ) + return mc_s + mc_a + + def compute_temperature_after_heating( + self, + initial_conductor_temperature: Celsius, + heating_duration: Duration, + current: Ampere, + time_step: Duration = np.timedelta64(1, "s"), + ) -> Celsius: + r"""Compute conductor temperature after it is heated for a time using an iterative method. + + Parameters + ---------- + initial_conductor_temperature: + :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. + heating_duration: + :math:`t`. Duration of heating. Unit is changed to seconds in the function. + current: + :math:`I~\left[\text{A}\right]`. Current heating the conductor + time_step: + :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`T~\left[^\circ\text{C}\right]`. Conductor temperature. + """ + def conductor_heating(temperature): + heat_capacity_ = self.compute_heat_capacity_per_unit_length(temperature) + heat_balance = self.compute_heat_balance(temperature, current=current) + return heat_balance / heat_capacity_ + # Convert durations to floating point numbers + dt = time_step / np.timedelta64(1, "s") + tf = heating_duration / np.timedelta64(1, "s") + final_temperature = solve_ivp_forward_euler(conductor_heating, initial_conductor_temperature, tf, dt) + return final_temperature + + def compute_transient_ampacity( + self, + max_conductor_temperature: Celsius, + heating_duration: Duration, + initial_conductor_temperature: Celsius, + time_step: Duration = np.timedelta64(60, "s"), + min_ampacity: Ampere = 0, + max_ampacity: Ampere = 5000, + tolerance: float = 1.0, + accept_invalid_values: bool = False, + ): + r"""Use the bisection method to compute the transient thermal rating (ampacity). + + Parameters + ---------- + max_conductor_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature + heating_duration: + :math:`t`. Duration of heating. Unit is changed to seconds in the function. + initial_conductor_temperature: + :math:`T_\text{init}~\left[^\circ\text{C}\right]`. Temperature of the conductor at the start of heating. + time_step: + :math:`\Delta t`. Time step of heating. Unit is changed to seconds in the function. + min_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for + computing the ampacity + max_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for + computing the ampacity + tolerance: + :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The + bisection iterations will stop once the numerical ampacity uncertainty is below + :math:`\Delta I`. The bisection method will run for + :math:`\left\lceil\frac{I_\text{min} - I_\text{min}}{\Delta I}\right\rceil` iterations. + accept_invalid_values: + If True, np.nan is returned whenever the current cannot be found within the provided + search interval. If False, a ValueError will be raised instead. + + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The transient thermal rating. + """ + n = self.span.num_conductors + I = solver.compute_conductor_transient_ampacity( # noqa + partial(self.compute_temperature_after_heating, time_step=time_step), + max_conductor_temperature, + initial_conductor_temperature, + heating_duration, + min_ampacity, + max_ampacity, + tolerance, + accept_invalid_values, + ) + return I * n diff --git a/linerate/solver.py b/linerate/solver.py index e6fc755..164850e 100644 --- a/linerate/solver.py +++ b/linerate/solver.py @@ -3,7 +3,7 @@ import numpy as np -from .units import Ampere, Celsius, FloatOrFloatArray, WattPerMeter +from .units import Ampere, Celsius, Duration, FloatOrFloatArray, WattPerMeter __all__ = ["bisect", "compute_conductor_temperature", "compute_conductor_ampacity"] @@ -164,3 +164,95 @@ def compute_conductor_ampacity( return bisect( f, min_ampacity, max_ampacity, tolerance, accept_invalid_values=accept_invalid_values ) + + +def compute_conductor_transient_ampacity( + final_temperature: Callable[[Celsius, Duration, Ampere], WattPerMeter], + max_conductor_temperature: Celsius, + initial_conductor_temperature: Celsius, + heating_duration: Duration, + min_ampacity: Ampere = 0, + max_ampacity: Ampere = 5_000, + tolerance: float = 1, # Ampere + accept_invalid_values: bool = False, +) -> Ampere: + r"""Use the bisection method to compute the temporary thermal rating (ampacity). + + Parameters + ---------- + final_temperature: + :math:`f(T, A) = P_J + P_s - P_c - P_r~\left[\text{W}~\text{m}^{-1}\right]`. A function of + temperature, heating time and current that returns the conductor temperature at the end of heating time. + initial_conductor_temperature: + :math:`T_\text{inital}~\left[^\circ\text{C}\right]`. Initial conductor temperature + heating_duration: + :math:`\delta t_\text{heating}~\left[\text{s}\right]`. Duration the conductor is heated. + max_conductor_temperature: + :math:`T_\text{max}~\left[^\circ\text{C}\right]`. Maximum allowed conductor temperature + min_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Lower bound for the numerical scheme for + computing the ampacity + max_ampacity: + :math:`I_\text{min}~\left[\text{A}\right]`. Upper bound for the numerical scheme for + computing the ampacity + tolerance: + :math:`\Delta I~\left[\text{A}\right]`. The numerical accuracy of the ampacity. The + bisection iterations will stop once the numerical ampacity uncertainty is below + :math:`\Delta I`. The bisection method will run for + :math:`\left\lceil\frac{I_\text{max} - I_\text{min}}{\Delta I}\right\rceil` iterations. + accept_invalid_values: + If True, np.nan is returned whenever the current cannot be found within the provided + search interval. If False, a ValueError will be raised instead. + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + :math:`I~\left[\text{A}\right]`. The thermal rating. + """ + + def temperature_difference(current: Ampere) -> Celsius: + return max_conductor_temperature - final_temperature( + initial_conductor_temperature, heating_duration, current + ) + + return bisect( + temperature_difference, + min_ampacity, + max_ampacity, + tolerance, + accept_invalid_values=accept_invalid_values, + ) + + +def solve_ivp_forward_euler( + f: Callable[[FloatOrFloatArray], FloatOrFloatArray], + y0: FloatOrFloatArray, + duration: FloatOrFloatArray, + step: FloatOrFloatArray, + ) -> FloatOrFloatArray: + """ + A very basic vectorized forward Euler solver time-invariant systems. + Parameters + ---------- + f: Time derivative of the system, assumed to only depend on state. dy / dt = f(y) + y0: Initial state of the system + duration: Time for which to apply the function + step: Time step + + Returns + ------- + Union[float, float64, ndarray[Any, dtype[float64]]] + Approximation of function f at the end of the given duration. + """ + y = y0 + step = np.broadcast_to(step, np.shape(duration)) + step_count = duration // step + remainder = duration % step + modification_mask = step_count > 0 + while np.any(modification_mask): + y = y + f(y) * step + step_count -= 1 + modification_mask = step_count > 0 + if np.any(remainder > 0): + y = y + f(y) * remainder + return y diff --git a/linerate/types.py b/linerate/types.py index 599e158..bbfc1ed 100644 --- a/linerate/types.py +++ b/linerate/types.py @@ -8,18 +8,29 @@ from .units import ( Celsius, Degrees, + JoulePerKilogramPerKelvin, + KilogramPerMeter, Meter, MeterPerSecond, OhmPerMeter, + PerKelvin, Radian, SquareMeter, SquareMeterPerAmpere, Unitless, - WattPerMeterPerKelvin, WattPerSquareMeter, + WattPerMeterPerKelvin, ) -__all__ = ["Conductor", "Weather", "Tower", "Span", "WeatherWithSolarRadiation"] +__all__ = [ + "Conductor", + "BaseWeather", + "Weather", + "Tower", + "Span", + "WeatherWithSolarRadiation", + "ConductorWithTransientData", +] @dataclass(frozen=True) @@ -75,6 +86,27 @@ class Conductor: thermal_conductivity: Optional[WattPerMeterPerKelvin] = None +@dataclass(frozen=True, kw_only=True) +class ConductorWithTransientData(Conductor): + #: :math:`m_s`. Mass of steel (core material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]` + steel_mass_per_unit_length: KilogramPerMeter + + #: :math:`c_s^20`. Heat capacity of steel (core material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]` + steel_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + + #: :math:`\beta_20^s`. Heat capacity temperature coefficient of the steel (core material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]` + steel_specific_heat_capacity_temperature_coefficient: PerKelvin + + #: :math:`m_a`. Mass of aluminium (outer material) per unit length of conductor :math:`\left[\text{kg}\text{m}^{-1}\right]` + aluminum_mass_per_unit_length: KilogramPerMeter + + #: :math:`c_a^20`. Heat capacity of aluminium (outer material) at 20 ^\circ C :math:`\left[\text{J}\text{kg}^{-1}\text{K}^{-1}\right]` + aluminum_specific_heat_capacity_at_20_celsius: JoulePerKilogramPerKelvin + + #: :math:`\beta_a^20`. Heat capacity temperature coefficient of the aluminium (outer material) at 20 ^\circ C :math:`\left[text{K}^{-1}\right]` + aluminum_specific_heat_capacity_temperature_coefficient: PerKelvin + + @dataclass(frozen=True) class Tower: """Container for a tower (span end-point).""" diff --git a/linerate/units.py b/linerate/units.py index af44c46..a833aae 100644 --- a/linerate/units.py +++ b/linerate/units.py @@ -32,5 +32,9 @@ SquareMeterPerSecond = Annotated[FloatOrFloatArray, "m²/s"] KilogramPerMeterPerSecond = Annotated[FloatOrFloatArray, "kg/(m s)"] KilogramPerCubeMeter = Annotated[FloatOrFloatArray, "kg/m³"] +KilogramPerMeter = Annotated[FloatOrFloatArray, "kg/m"] +PerKelvin = Annotated[FloatOrFloatArray, "1/K"] +JoulePerMeter = Annotated[FloatOrFloatArray, "J/m"] Date = Union[np.datetime64, npt.NDArray[np.datetime64]] +Duration = Union[np.timedelta64, npt.NDArray[np.timedelta64]] diff --git a/pyproject.toml b/pyproject.toml index 4e0368c..aa8035a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -84,7 +84,7 @@ source = ["linerate"] log_file = "pytest.log" log_level = "DEBUG" log_file_format = "%(asctime)s - %(name)s - %(levelname)s: %(message)s" -norecursedirs = "_build tmp* __pycache__ src prof wheel_files" +norecursedirs = ".hypothesis _build tmp* __pycache__ src prof wheel_files" markers = [ "integration: Mark test as an integration test", ] diff --git a/tests/acceptance_tests/test_cigre_transient_cases.py b/tests/acceptance_tests/test_cigre_transient_cases.py new file mode 100644 index 0000000..6ec6e7c --- /dev/null +++ b/tests/acceptance_tests/test_cigre_transient_cases.py @@ -0,0 +1,119 @@ +import numpy as np +import pytest +from pytest import approx +from linerate import ConductorWithTransientData, Weather, Tower, Span, Cigre601 + + +@pytest.fixture +def drake_conductor_transient() -> ConductorWithTransientData: + return ConductorWithTransientData( + conductor_diameter=28.143e-3, + core_diameter=10.4e-3, + outer_layer_strand_diameter=4.4e-3, + emissivity=0.8, + solar_absorptivity=0.8, + temperature1=25, + temperature2=75, + resistance_at_temperature1=0.0727e-3, + resistance_at_temperature2=0.0872e-3, + aluminium_cross_section_area=float("nan"), + constant_magnetic_effect=1, + current_density_proportional_magnetic_effect=0, + max_magnetic_core_relative_resistance_increase=1, + steel_mass_per_unit_length=0.5119, + steel_specific_heat_capacity_at_20_celsius=481, + steel_specific_heat_capacity_temperature_coefficient=1e-4, + aluminum_mass_per_unit_length=1.116, + aluminum_specific_heat_capacity_at_20_celsius=897, + aluminum_specific_heat_capacity_temperature_coefficient=3.8e-4, + ) + + +@pytest.fixture +def initial_weather() -> Weather: + return Weather( + air_temperature=24, + wind_speed=1.9, + wind_direction=np.radians(35), + ground_albedo=1.0, + ) + + +@pytest.fixture +def weather_1() -> Weather: + return Weather( + air_temperature=23.7, + wind_speed=1.7, + wind_direction=np.radians(28), + ground_albedo=1.0, + ) + + +@pytest.fixture +def weather_2() -> Weather: + return Weather( + air_temperature=23.5, + wind_speed=0.8, + wind_direction=np.radians(53), + ground_albedo=1.0, + ) + + +@pytest.fixture +def example_span_transient(drake_conductor_transient: ConductorWithTransientData) -> Span: + start_tower = Tower(latitude=50, longitude=0.0001, altitude=0) + end_tower = Tower(latitude=50, longitude=-0.0001, altitude=0) + return Span( + conductor=drake_conductor_transient, + start_tower=start_tower, + end_tower=end_tower, + num_conductors=1, + ) + + +@pytest.fixture +def transient_initial_model(example_span_transient, initial_weather) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=initial_weather, time=np.datetime64("2016-06-10 00:00") + ) + + +@pytest.fixture +def transient_model_1(example_span_transient, weather_1) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=weather_1, time=np.datetime64("2016-06-10 00:10") + ) + + +@pytest.fixture +def transient_model_2(example_span_transient, weather_2) -> Cigre601: + return Cigre601( + span=example_span_transient, weather=weather_2, time=np.datetime64("2016-06-10 00:20") + ) + + +def test_initial_condition(transient_initial_model): + assert transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) == approx( + 42.01, abs=0.01 + ) + + +def test_tracking(transient_initial_model, transient_model_1, transient_model_2): + start_temperature = transient_initial_model.compute_conductor_temperature(802, tolerance=0.01) + # 60s time step is used since the calculations in the brochure use the same step size. + results_1 = [ + transient_model_1.compute_temperature_after_heating( + start_temperature, np.timedelta64(t, "m"), 819, time_step=np.timedelta64(60, "s") + ) + for t in range(1, 11) + ] + results_2 = [ + transient_model_2.compute_temperature_after_heating( + results_1[-1], np.timedelta64(t, "m"), 856, time_step=np.timedelta64(60, "s") + ) + for t in range(1, 11) + ] + correct_1 = [42.175, 42.321, 42.449, 42.562, 42.662, 42.750, 42.828, 42.897, 42.958, 43.011] + correct_2 = [44.147, 45.199, 46.174, 47.075, 47.910, 48.682, 49.396, 50.057, 50.668, 51.233] + assert all(result == approx(correct, abs=0.01) for result, correct in zip(results_1, correct_1)) + assert all(result == approx(correct, abs=0.01) for result, correct in zip(results_2, correct_2)) diff --git a/tests/equations/test_heat_capacity.py b/tests/equations/test_heat_capacity.py new file mode 100644 index 0000000..fe1c2ea --- /dev/null +++ b/tests/equations/test_heat_capacity.py @@ -0,0 +1,55 @@ +from pytest import approx + +import hypothesis +import hypothesis.strategies as st + +from linerate.units import Celsius, KilogramPerMeter +from linerate.equations import heat_capacity + + +@hypothesis.given(mass_per_unit_length=st.floats(min_value=0, max_value=1e10, allow_nan=False)) +def test_heat_capacity_is_linear_in_mass_per_unit_length(mass_per_unit_length: KilogramPerMeter): + m = mass_per_unit_length + c = 1 + beta = 1 + T = 20 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(m) + + +@hypothesis.given( + specific_heat_capacity_at_20_celsius=st.floats(min_value=0, max_value=1e10, allow_nan=False) +) +def test_heat_capacity_is_linear_in_specific_heat_capacity( + specific_heat_capacity_at_20_celsius: KilogramPerMeter, +): + m = 1 + c = specific_heat_capacity_at_20_celsius + beta = 1 + T = 20 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(c) + + +@hypothesis.given( + specific_heat_capacity_coefficient=st.floats(min_value=0, max_value=1e10, allow_nan=False) +) +def test_heat_capacity_is_linear_in_heat_capacity_coefficient( + specific_heat_capacity_coefficient: KilogramPerMeter, +): + m = 1 + c = 1 + beta = specific_heat_capacity_coefficient + T = 21 + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(beta + 1) + + +@hypothesis.given(temperature=st.floats(min_value=-1e10, max_value=1e10, allow_nan=False)) +def test_heat_capacity_is_linear_in_temperature(temperature: Celsius): + m = 1 + c = 1 + beta = 1 + T = temperature + mc = heat_capacity.calculate_heat_capacity_per_unit_length(m, c, beta, T) + assert mc == approx(T - 19) diff --git a/tests/integration_tests/test_vectorization.py b/tests/integration_tests/test_vectorization.py index 4adb809..2ac2ee1 100644 --- a/tests/integration_tests/test_vectorization.py +++ b/tests/integration_tests/test_vectorization.py @@ -1,67 +1,127 @@ import hypothesis -import hypothesis.strategies as st +from hypothesis.extra.numpy import mutually_broadcastable_shapes, BroadcastableShapes import numpy as np +import numpy.typing as npt from linerate.models.cigre601 import Cigre601 -from linerate.types import Conductor, Span, Tower, Weather +from linerate.types import Conductor, ConductorWithTransientData, Span, Tower, Weather + +NUMBER_OF_VARIABLES_TRANSIENT = 33 +NUMBER_OF_VARIABLES_STEADY_STATE = 24 + + +class Shapes: + def __init__(self, shapes: tuple[tuple[int, ...], ...]): + self.shapes = shapes + self.index = 0 + + def get_zeros(self, dtype: np.dtype = np.dtype("f8")) -> npt.NDArray: + out = np.zeros(self.shapes[self.index], dtype=dtype) + self.index += 1 + return out @hypothesis.given( - vectorization_indices=st.sets(st.integers(min_value=1, max_value=24), max_size=10) + array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES_STEADY_STATE, max_dims=5) ) -def test_vectorization(vectorization_indices: set[int]): - def get_zeros(index, length): - out = np.ones(length, dtype=int) - if index in vectorization_indices: - out[index] = 2 - return np.zeros(out) +def test_vectorization(array_shapes: BroadcastableShapes): + shapes = Shapes(array_shapes.input_shapes) conductor = Conductor( - core_diameter=get_zeros(0, 25) + 1, - conductor_diameter=get_zeros(1, 25) + 2, - outer_layer_strand_diameter=get_zeros(2, 25), - emissivity=get_zeros(3, 25), - solar_absorptivity=get_zeros(4, 25), - temperature1=get_zeros(5, 25), - temperature2=get_zeros(6, 25) + 1, - resistance_at_temperature1=get_zeros(7, 25), - resistance_at_temperature2=get_zeros(8, 25), - aluminium_cross_section_area=1 + get_zeros(9, 25), - constant_magnetic_effect=get_zeros(10, 25), - current_density_proportional_magnetic_effect=get_zeros(11, 25) + 1, - max_magnetic_core_relative_resistance_increase=get_zeros(12, 25), + core_diameter=1, + conductor_diameter=shapes.get_zeros() + 2, + outer_layer_strand_diameter=shapes.get_zeros(), + emissivity=shapes.get_zeros(), + solar_absorptivity=shapes.get_zeros(), + temperature1=shapes.get_zeros(), + temperature2=shapes.get_zeros() + 1, + resistance_at_temperature1=shapes.get_zeros(), + resistance_at_temperature2=shapes.get_zeros(), + aluminium_cross_section_area=1 + shapes.get_zeros(), + constant_magnetic_effect=shapes.get_zeros(), + current_density_proportional_magnetic_effect=shapes.get_zeros() + 1, + max_magnetic_core_relative_resistance_increase=shapes.get_zeros(), ) weather = Weather( - air_temperature=get_zeros(13, 25), - wind_direction=get_zeros(14, 25), - wind_speed=get_zeros(15, 25), - ground_albedo=get_zeros(23, 25), - clearness_ratio=get_zeros(16, 25), + air_temperature=shapes.get_zeros(), + wind_direction=shapes.get_zeros(), + wind_speed=shapes.get_zeros(), + ground_albedo=shapes.get_zeros(), + clearness_ratio=shapes.get_zeros(), ) span = Span( conductor=conductor, start_tower=Tower( - get_zeros(17, 25), - get_zeros(18, 25), - get_zeros(19, 25), + shapes.get_zeros(), + shapes.get_zeros(), + shapes.get_zeros(), ), end_tower=Tower( - 0.1 + get_zeros(20, 25), - 0.1 + get_zeros(21, 25), - 10 + get_zeros(22, 25), + 0.1 + shapes.get_zeros(), + 0.1 + shapes.get_zeros(), + 10 + shapes.get_zeros(), ), num_conductors=1, ) - time = np.array([np.datetime64("2022-01-01")]).reshape([1] * 25) - if 24 in vectorization_indices: - time = np.concatenate([time, time], axis=-1) + time = shapes.get_zeros(dtype=np.dtype("datetime64[D]")) model = Cigre601(span=span, weather=weather, time=time) temperature = model.compute_conductor_temperature(100) - shape = [1] * 25 - for vectorisation_index in vectorization_indices: - shape[vectorisation_index] = 2 - shape[0] = 1 - assert isinstance(temperature, np.ndarray) - assert temperature.shape == tuple(shape) + assert temperature.shape == array_shapes.result_shape + + +@hypothesis.given(array_shapes=mutually_broadcastable_shapes(num_shapes=NUMBER_OF_VARIABLES_TRANSIENT, max_dims=5)) +def test_transient_vectorization(array_shapes: BroadcastableShapes): + shapes = Shapes(array_shapes.input_shapes) + + conductor = ConductorWithTransientData( + core_diameter=1, # Core diameter is not used in any calculations so its shape will not be broadcast + conductor_diameter=shapes.get_zeros() + 2, + outer_layer_strand_diameter=shapes.get_zeros(), + emissivity=shapes.get_zeros(), + solar_absorptivity=shapes.get_zeros(), + temperature1=shapes.get_zeros(), + temperature2=shapes.get_zeros() + 1, + resistance_at_temperature1=shapes.get_zeros(), + resistance_at_temperature2=shapes.get_zeros(), + aluminium_cross_section_area=1 + shapes.get_zeros(), + constant_magnetic_effect=shapes.get_zeros(), + current_density_proportional_magnetic_effect=shapes.get_zeros() + 1, + max_magnetic_core_relative_resistance_increase=shapes.get_zeros(), + steel_mass_per_unit_length=shapes.get_zeros() + 0.5, + steel_specific_heat_capacity_at_20_celsius=shapes.get_zeros() + 500, + steel_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), + aluminum_mass_per_unit_length=shapes.get_zeros() + 1, + aluminum_specific_heat_capacity_at_20_celsius=shapes.get_zeros() + 1000, + aluminum_specific_heat_capacity_temperature_coefficient=shapes.get_zeros(), + ) + weather = Weather( + air_temperature=shapes.get_zeros(), + wind_direction=shapes.get_zeros(), + wind_speed=shapes.get_zeros(), + ground_albedo=shapes.get_zeros(), + clearness_ratio=shapes.get_zeros(), + ) + span = Span( + conductor=conductor, + start_tower=Tower( + shapes.get_zeros(), + shapes.get_zeros(), + shapes.get_zeros(), + ), + end_tower=Tower( + 0.1 + shapes.get_zeros(), + 0.1 + shapes.get_zeros(), + 10 + shapes.get_zeros(), + ), + num_conductors=1, + ) + time = shapes.get_zeros(dtype=np.dtype("datetime64[D]")) + model = Cigre601(span=span, weather=weather, time=time) + final_temperature = model.compute_temperature_after_heating( + shapes.get_zeros() + 20, + shapes.get_zeros(dtype=np.dtype("