From f243ced992fbdfe8243323cdc92d0ee97abf967f Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:02:38 +0200 Subject: [PATCH 01/53] Changed default value of srad parameter in all solar heating power to nan instead of None; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/power/cigre/solar_heating.py | 9 +++++---- src/thermohl/power/cner/solar_heating.py | 7 ++++--- src/thermohl/power/ieee/solar_heating.py | 7 ++++--- src/thermohl/power/olla/solar_heating.py | 9 +++++---- test/test_power_cner_integration.py | 2 +- 5 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/thermohl/power/cigre/solar_heating.py b/src/thermohl/power/cigre/solar_heating.py index d4482a8..eeed150 100644 --- a/src/thermohl/power/cigre/solar_heating.py +++ b/src/thermohl/power/cigre/solar_heating.py @@ -48,7 +48,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -74,8 +74,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx. + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns @@ -85,7 +86,7 @@ def __init__( """ self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = SolarHeating._solar_radiation( np.deg2rad(lat), np.deg2rad(azm), al, month, day, hour ) diff --git a/src/thermohl/power/cner/solar_heating.py b/src/thermohl/power/cner/solar_heating.py index 74b6e89..0e8a053 100644 --- a/src/thermohl/power/cner/solar_heating.py +++ b/src/thermohl/power/cner/solar_heating.py @@ -21,7 +21,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Build with args. @@ -45,8 +45,9 @@ def __init__( external diameter. alpha : np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/ieee/solar_heating.py b/src/thermohl/power/ieee/solar_heating.py index 67008c5..a1c66c8 100644 --- a/src/thermohl/power/ieee/solar_heating.py +++ b/src/thermohl/power/ieee/solar_heating.py @@ -25,7 +25,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -53,8 +53,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/olla/solar_heating.py b/src/thermohl/power/olla/solar_heating.py index 1df8d69..99bc5a7 100644 --- a/src/thermohl/power/olla/solar_heating.py +++ b/src/thermohl/power/olla/solar_heating.py @@ -24,7 +24,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -38,7 +38,7 @@ def __init__( lat : float or np.ndarray Latitude. alt : float or np.ndarray - Lltitude. + Longitude. azm : float or np.ndarray Azimuth. month : int or np.ndarray @@ -52,8 +52,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index d41ca0f..58b7ea6 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -220,6 +220,6 @@ def test_solar_heating(): alpha = 0.9 * ones p = np.array([34.9, 21.9357, 13.95, 21.9357, 21.9357]) - s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha, srad=None) + s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha) assert np.allclose(p, s.value(ones), 0.1) From fa5730e0e98b32a1abfa14f2c629c83a0f0dc021 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:05:38 +0200 Subject: [PATCH 02/53] Changed default value of srad parameter in all base solar heating power term to nan instead of None Signed-off-by: Emmanuel Cieren --- src/thermohl/power/solar_heating.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/thermohl/power/solar_heating.py b/src/thermohl/power/solar_heating.py index f4a9c90..72ec76b 100644 --- a/src/thermohl/power/solar_heating.py +++ b/src/thermohl/power/solar_heating.py @@ -66,11 +66,11 @@ def __init__( D: floatArrayLike, alpha: floatArrayLike, est: _SRad, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = est(np.deg2rad(lat), alt, np.deg2rad(azm), tb, month, day, hour) else: self.srad = np.maximum(srad, 0.0) From b79c97601743118074859a03de8cb792e7a162b5 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:06:25 +0200 Subject: [PATCH 03/53] Changed some var names to more readable, private values Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/__init__.py | 40 ++++++++++++++++----------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index b987948..36d2f30 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -9,10 +9,10 @@ from typing import Dict, Any, Optional, Union, Type -from thermohl.power import cigre as cigrep -from thermohl.power import cner as cnerp -from thermohl.power import ieee as ieeep -from thermohl.power import olla as ollap +from thermohl.power import cigre as _cigre +from thermohl.power import cner as _cner +from thermohl.power import ieee as _ieee +from thermohl.power import olla as _olla from thermohl.solver.base import Args, Solver from thermohl.solver.slv1d import Solver1D @@ -45,34 +45,34 @@ def _factory( if model == "cigre": return solver( dic, - cigrep.JouleHeating, - cigrep.SolarHeating, - cigrep.ConvectiveCooling, - cigrep.RadiativeCooling, + _cigre.JouleHeating, + _cigre.SolarHeating, + _cigre.ConvectiveCooling, + _cigre.RadiativeCooling, ) elif model == "ieee": return solver( dic, - ieeep.JouleHeating, - ieeep.SolarHeating, - ieeep.ConvectiveCooling, - ieeep.RadiativeCooling, + _ieee.JouleHeating, + _ieee.SolarHeating, + _ieee.ConvectiveCooling, + _ieee.RadiativeCooling, ) elif model == "olla": return solver( dic, - ollap.JouleHeating, - ollap.SolarHeating, - ollap.ConvectiveCooling, - ollap.RadiativeCooling, + _olla.JouleHeating, + _olla.SolarHeating, + _olla.ConvectiveCooling, + _olla.RadiativeCooling, ) elif model == "cner": return solver( dic, - cnerp.JouleHeating, - cnerp.SolarHeating, - cnerp.ConvectiveCooling, - cnerp.RadiativeCooling, + _cner.JouleHeating, + _cner.SolarHeating, + _cner.ConvectiveCooling, + _cner.RadiativeCooling, ) else: raise ValueError() From d960b77623a4da42a3cf186e51b079a378221b1c Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:09:14 +0200 Subject: [PATCH 04/53] Added solar irradiance (srad) to Solver.args; updated comments Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index bc13416..1424a1f 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -62,9 +62,9 @@ def _set_default_values(self) -> None: self.alt = 0.0 # altitude (m) self.azm = 0.0 # azimuth (deg) - self.month = 3 # month number (1=Jan, 2=Feb, ...) + self.month = 3 # month number (1=Jan, 2=Feb ...) self.day = 21 # day of the month - self.hour = 12 # hour of the day (in [0, 23] range) + self.hour = 12.0 # hour of the day (in [0, 23] range) self.Ta = 15.0 # ambient temperature (C) self.Pa = 1.0e05 # ambient pressure (Pa) @@ -73,8 +73,8 @@ def _set_default_values(self) -> None: self.ws = 0.0 # wind speed (m.s**-1) self.wa = 90.0 # wind angle (deg, regarding north) self.al = 0.8 # albedo (1) - # coefficient for air pollution from 0 (clean) to 1 (polluted) - self.tb = 0.1 + self.tb = 0.1 # turbidity (coefficient for air pollution from 0 -clean- to 1 -polluted-) + self.srad = float("nan") # solar irradiance (in W.m**-2) self.I = 100.0 # transit intensity (A) From 2ba2a50d1ee17aa9b7aa6f42296fd338df6e6e9f Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:11:18 +0200 Subject: [PATCH 05/53] Renamed function with more explicit name; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/utils.py | 2 +- test/test_power_cner_integration.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index af82538..5fa21ac 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -101,7 +101,7 @@ def add_default_uncertainties(dat: dict, warning: bool = False) -> dict: return _dict_completion(dat, fil, check=False, warning=warning) -def df2dct(df: pd.DataFrame) -> dict: +def df2dict(df: pd.DataFrame) -> dict: """Convert a pandas.DataFrame to a dictionary. Would be an equivalent to df.to_dict(orient='numpy.ndarray') if it existed. diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index 58b7ea6..3f924cc 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -184,13 +184,13 @@ def test_compare_power(): ds = pd.concat(len(T) * (ds,)).reset_index(drop=True) T = np.concatenate([n * (t,) for t in T]) - from thermohl.utils import df2dct + from thermohl.utils import df2dict - d1 = df2dct(ds) + d1 = df2dict(ds) ds["wa"] = np.rad2deg( np.arcsin(np.sin(np.deg2rad(np.abs(ds["azm"] - ds["wa"]) % 180.0))) ) - d2 = df2dct(ds) + d2 = df2dict(ds) del (ds, n) pj = cner.JouleHeating(**d1) From 86ab48a9345cd882a163dc1a319d512018254516 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:53:12 +0200 Subject: [PATCH 06/53] Changed max_len/exten/compress methods in Solver class; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 98 +++++++++++++++-------------- src/thermohl/solver/slv1t.py | 11 +++- src/thermohl/solver/slv3t.py | 14 ++--- src/thermohl/solver/slv3t_legacy.py | 10 +-- test/test_shape.py | 2 +- test/unit/solver/test_base.py | 56 ++++++++--------- 6 files changed, 101 insertions(+), 90 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 1424a1f..4ec09f5 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -9,7 +9,7 @@ import datetime from abc import ABC, abstractmethod -from typing import Tuple, Type, Any, Optional, KeysView, Dict +from typing import Tuple, Type, Any, Optional, Dict, KeysView import numpy as np import pandas as pd @@ -53,6 +53,8 @@ def __init__(self, dic: Optional[dict[str, Any]] = None): for k in dic: if k in keys and dic[k] is not None: self[k] = dic[k] + # check for shape incompatibilities + _ = self.shape() def _set_default_values(self) -> None: """Set default values.""" @@ -115,53 +117,52 @@ def __getitem__(self, key: str) -> Any: def __setitem__(self, key: str, value: Any) -> None: self.__dict__[key] = value - def max_len(self) -> int: + def shape(self) -> Tuple[int, ...]: """ - Calculate the maximum length of the values in the dictionary. + Compute the maximum effective shape of values in current instance. - This method iterates over all keys in the dictionary and determines the maximum length - of the values associated with those keys. - If a value is not of a type that has a length, it is ignored. + Members of Args can be float of arrays. If arrays, they must be + one-diemnsional. Float and 1d array can coexist, but all arrays should + have the same shape/size. - Returns: - int: The maximum length of the values in the dictionary. If the dictionary is empty - or all values are of types that do not have a length, the method returns 1. - """ - n = 1 - for k in self.keys(): - try: - n = max(n, len(self[k])) - except TypeError: - pass - return n - - def extend_to_max_len(self) -> None: - """ - Extend all elements in the dictionary to the maximum length. + This method iterates over all keys in the instance's __dict__ and + computes the maximum length of the values associated with those keys. - This method iterates over all keys in the dictionary and checks if the - corresponding value is a numpy ndarray. If it is, it checks if its length - matches the maximum length obtained from the `max_len` method. - If the length matches, it creates a copy of the array. - If the length does not match or for non-ndarray values, it creates - a new numpy array of the maximum length, filled with the original value - and having the same data type. + If incompatible shapes are encountered, an exception is raised (ValueError). - Returns: - None """ - n = self.max_len() + shape_ = () for k in self.keys(): - if isinstance(self[k], np.ndarray): - t = self[k].dtype - c = len(self[k]) == n - else: - t = type(self[k]) - c = False - if c: - self[k] = self[k][:] + s = np.array(self[k]).shape + d = len(s) + er = f"Key {k} has a {s} shape when main shape is {shape_}" + if d == 0: + continue + if d == 1: + if shape_ == (): + shape_ = s + elif len(shape_) == 1: + if shape_ != s: + raise ValueError(er) + else: + raise ValueError(er) else: - self[k] = self[k] * np.ones((n,), dtype=t) + raise ValueError( + f"Key {k} has a {s} shape, only float and 1-dim arrays are accepted" + ) + return shape_ + + def extend(self, shape: Tuple[int, ...] = None) -> None: + # get shape + if shape is None: + shape = self.shape() + # complete if necessary + if len(shape) == 0: + shape = (1,) + # create a copy dict with scaled array + for k in self.keys(): + a = np.array(self[k]) + self[k] = a * np.ones(shape, dtype=a.dtype) def compress(self) -> None: """ @@ -172,10 +173,9 @@ def compress(self) -> None: None """ for k in self.keys(): - if isinstance(self[k], np.ndarray): - u = np.unique(self[k]) - if len(u) == 1: - self[k] = u[0] + u = np.unique(self[k]) + if len(u) == 1: + self[k] = u[0] class Solver(ABC): @@ -246,7 +246,7 @@ def __init__( """ self.args = Args(dic) - self.args.extend_to_max_len() + self.args.extend() self.jh = joule(**self.args.__dict__) self.sh = solar(**self.args.__dict__) self.cc = convective(**self.args.__dict__) @@ -254,8 +254,14 @@ def __init__( self.pc = precipitation(**self.args.__dict__) self.args.compress() + def _min_shape(self) -> Tuple[int, ...]: + shape = self.args.shape() + if shape == (): + shape = (1,) + return shape + def update(self) -> None: - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 04f735f..39f9bcc 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -57,7 +57,12 @@ def steady_temperature( # solve with bisection T, err = bisect_v( - lambda x: -self.balance(x), Tmin, Tmax, (self.args.max_len(),), tol, maxiter + lambda x: -self.balance(x), + Tmin, + Tmax, + self._min_shape(), + tol=tol, + maxiter=maxiter, ) # format output @@ -102,7 +107,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError("The length of the time array must be at least 2.") @@ -223,7 +228,7 @@ def steady_intensity( transit = self.args.I # solve with bisection - shape = (self.args.max_len(),) + shape = self._min_shape() T_ = T * np.ones(shape) jh = ( self.cc.value(T_) diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index f038a0b..d04b790 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -78,7 +78,7 @@ def _morgan_coefficients( - i : numpy.ndarray[int] Indices where surface diameter `d_` is greater than 0. """ - c = 0.5 * np.ones((self.args.max_len(),)) + c = 0.5 * np.ones(self._min_shape()) D = self.args.D * np.ones_like(c) d = self.args.d * np.ones_like(c) i = np.nonzero(d > 0.0)[0] @@ -97,7 +97,7 @@ def update(self) -> None: Returns: None """ - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) @@ -215,7 +215,7 @@ def steady_temperature( """ # if no guess provided, use ambient temp - shape = (self.args.max_len(),) + shape = self._min_shape() Tsg = Tsg if Tsg is not None else 1.0 * self.args.Ta Tcg = Tcg if Tcg is not None else 1.5 * np.abs(self.args.Ta) Tsg_ = Tsg * np.ones(shape) @@ -322,7 +322,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() @@ -436,9 +436,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes c, D, d, ix = self.mgc diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 2da3a2f..561dfcb 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -46,7 +46,7 @@ def _morgan_coefficients(self) -> Tuple[floatArray, intArray]: UNIFORM_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 13 BIMETALLIC_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 21 - core_diameter_array = self.args.d * np.ones((self.args.max_len(),)) + core_diameter_array = self.args.d * np.ones(self._min_shape()) indices_non_zero_diameter = np.nonzero(core_diameter_array > 0.0)[0] heat_flux_coefficients = UNIFORM_CONDUCTOR_COEFFICIENT * np.ones_like( core_diameter_array @@ -91,9 +91,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -162,7 +162,7 @@ def transient_temperature_legacy( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() diff --git a/test/test_shape.py b/test/test_shape.py index 0c1f1ab..20be2f4 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -95,7 +95,7 @@ def test_steady_1d_mix(): n = 61 for s in _solvers(): s.args.Ta = np.linspace(-10, +50, n) - s.args.I = np.array([199.0]) + s.args.I = np.array(199.0) s.update() t = s.steady_temperature() a = _ampargs(s, t) diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 63eb841..9aba5d1 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -16,59 +16,59 @@ def test_max_len_with_mixed_types(): dic = { "lat": np.array([45.0, 46.0]), "lon": 10.0, - "alt": np.array([20.0, 25.0, 30.0]), + "alt": np.array([20.0, 25.0]), } args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 3 + assert result == (2,) -def test_max_len_with_ndarray(): +def test_shape_with_ndarray(): dic = {"lat": np.array([45.0, 46.0]), "lon": np.array([10.0, 11.0])} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 2 + assert result == (2,) -def test_max_len_with_scalar(): +def test_shape_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_empty_dict(): +def test_shape_with_empty_dict(): args = Args({}) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_varied_lengths(): +def test_shape_with_varied_lengths(): dic = { "lat": np.array([45.0, 46.0]), "lon": np.array([10.0]), "alt": np.array([20.0, 25.0, 30.0]), } - args = Args(dic) - - result = args.max_len() + try: + args = Args(dic) + assert False + except ValueError: + pass - assert result == 3 - -def test_extend_to_max_len_with_nd_array(): +def test_extend_with_nd_array(): dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -78,11 +78,11 @@ def test_extend_to_max_len_with_nd_array(): np.testing.assert_array_equal(args.lon, np.array([10.0, 10.0])) -def test_extend_to_max_len_with_scalar(): +def test_extend_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -92,11 +92,11 @@ def test_extend_to_max_len_with_scalar(): np.testing.assert_array_equal(args.lon, np.array([10.0])) -def test_extend_to_max_len_with_mixed_types(): - dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array([20.0])} +def test_extend_with_mixed_types(): + dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array(20.0)} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -109,10 +109,10 @@ def test_extend_to_max_len_with_mixed_types(): np.testing.assert_array_equal(args.alt, np.array([20.0, 20.0])) -def test_extend_to_max_len_with_empty_dict(): +def test_extend_with_empty_dict(): args = Args({}) - args.extend_to_max_len() + args.extend() for key in args.keys(): assert isinstance(args[key], (float, int, ndarray)) @@ -164,7 +164,7 @@ def test_compress_with_empty_dict(): args.compress() for key in args.keys(): - assert isinstance(args[key], (float, int, ndarray)) + assert isinstance(args[key], (float, int, np.integer, np.floating, ndarray)) if isinstance(args[key], ndarray): assert len(args[key]) == 1 From b6700b8f4b67ba076553bfcb455bf919badcf9b2 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:55:28 +0200 Subject: [PATCH 07/53] Fixed typo/error in set_times method; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 5 +---- test/unit/solver/test_base.py | 10 +++++----- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 4ec09f5..b1fe600 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -368,10 +368,7 @@ def _set_dates( td = np.array( [datetime.timedelta()] - + [ - datetime.timedelta(seconds=float(time[i] - time[i - 1])) - for i in range(1, N) - ] + + [datetime.timedelta(seconds=float(time[i] - time[0])) for i in range(1, N)] ) for j in range(n): diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 9aba5d1..089ed0c 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -256,7 +256,7 @@ def test_set_dates_single_day(): assert hours[1, 0] == 1.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 def test_set_dates_multiple_days(): @@ -279,7 +279,7 @@ def test_set_dates_multiple_days(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 2 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_months(): @@ -302,7 +302,7 @@ def test_set_dates_multiple_months(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_inputs(): @@ -328,7 +328,7 @@ def test_set_dates_multiple_inputs(): assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 assert months[0, 1] == 2 assert days[0, 1] == 2 @@ -340,4 +340,4 @@ def test_set_dates_multiple_inputs(): assert months[2, 1] == 2 assert days[2, 1] == 2 - assert hours[2, 1] == 13.0 + assert hours[2, 1] == 14.0 From b14a53b6eab924233db77757812756fec1044eee Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:58:08 +0200 Subject: [PATCH 08/53] Changed reshape function to be more restrictive; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 39 ++++++++++++++++++++--------------- test/unit/solver/test_base.py | 15 ++++++++++++-- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index b1fe600..7034bfe 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -291,12 +291,12 @@ def steady_intensity(self) -> pd.DataFrame: raise NotImplementedError -def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: +def reshape(input_var: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: """ Reshape the input array to the specified dimensions (nr, nc) if possible. Args: - input_array (numberArrayLike): Input array to be reshaped. + input (numberArrayLike): Variable to be reshaped. nb_row (int): Desired number of rows for the reshaped array. nb_columns (int): Desired number of columns for the reshaped array. @@ -305,24 +305,29 @@ def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numbe returns an array filled with the input_value repeated to fill the dimension (nb_row, nb_columns). Raises: - AttributeError: If the input_array has an invalid shape that cannot be reshaped. + ValueError: If the input has an invalid shape that cannot be reshaped. """ - reshaped_array = ndarray - try: - input_shape = input_array.shape - if len(input_shape) == 1: - if nb_row == input_shape[0]: - reshaped_array = np.column_stack(nb_columns * (input_array,)) - elif nb_columns == input_shape[0]: - reshaped_array = np.vstack(nb_row * (input_array,)) - elif len(input_shape) == 0: - raise AttributeError() - else: - reshaped_array = np.reshape(input_array, (nb_row, nb_columns)) - except AttributeError: + + input_array = np.array(input_var) + input_shape = input_array.shape + + msg = f"Input array has incompatible shape {input_shape} with specified number of rows ({nb_row}) and/or columns ({nb_columns})." + + if len(input_shape) == 0: reshaped_array = input_array * np.ones( - (nb_row, nb_columns), dtype=type(input_array) + (nb_row, nb_columns), dtype=input_array.dtype ) + elif len(input_shape) == 1: + if nb_row == input_shape[0]: + reshaped_array = np.column_stack(nb_columns * (input_array,)) + elif nb_columns == input_shape[0]: + reshaped_array = np.vstack(nb_row * (input_array,)) + else: + raise ValueError(msg) + elif input_shape == (nb_row, nb_columns): + return input_array + else: + raise ValueError(msg) return reshaped_array diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 089ed0c..f3822aa 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -225,8 +225,8 @@ def test_reshape_scalar_to_2d(): np.testing.assert_array_equal(result, expected) -def test_reshape_invalid_shape(): - array = np.array(0) # ([1.0, 2.0, 3.0]) +def test_reshape_another_scalar_to_2d(): + array = np.array(0) nb_row = 2 nb_columns = 2 expected = np.array([[0, 0], [0, 0]]) @@ -236,6 +236,17 @@ def test_reshape_invalid_shape(): np.testing.assert_array_equal(result, expected) +def test_reshape_invalid_shape(): + array = np.array([1.0, 2.0, 3.0]) + nb_row = 2 + nb_columns = 2 + try: + reshape(array, nb_row, nb_columns) + assert False + except ValueError: + pass + + def test_set_dates_single_day(): month = np.array([1]) day = np.array([1]) From 442e3aeac6168df2570687f78603428e3beb4a77 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:07:46 +0200 Subject: [PATCH 09/53] Factorized some code for steady temperature methods. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 47 ++++++++++++++++++-------------- src/thermohl/solver/slv3t.py | 53 +++++++++++++++++++++--------------- 2 files changed, 58 insertions(+), 42 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 39f9bcc..a695a89 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -20,6 +20,31 @@ class Solver1T(Solver_): + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + T: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(T) + df[Solver_.Names.psol] = self.sh.value(T) + df[Solver_.Names.pcnv] = self.cc.value(T) + df[Solver_.Names.prad] = self.rc.value(T) + df[Solver_.Names.ppre] = self.pc.value(T) + + return df + def steady_temperature( self, Tmin: float = DP.tmin, @@ -67,16 +92,7 @@ def steady_temperature( # format output df = pd.DataFrame(data=T, columns=[Solver_.Names.temp]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T, err, df) return df @@ -249,15 +265,6 @@ def fun(i: floatArray) -> floatArrayLike: # format output df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T_, err, df) return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d04b790..d8bd60a 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -183,6 +183,32 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + Ts: np.ndarray, + Ta: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(Ta) + df[Solver_.Names.psol] = self.sh.value(Ts) + df[Solver_.Names.pcnv] = self.cc.value(Ts) + df[Solver_.Names.prad] = self.rc.value(Ts) + df[Solver_.Names.ppre] = self.pc.value(Ts) + + return df + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -240,16 +266,7 @@ def steady_temperature( df = pd.DataFrame( {Solver_.Names.tsurf: x, Solver_.Names.tavg: z, Solver_.Names.tcore: y} ) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.joule(x, y) - df[Solver_.Names.psol] = self.sh.value(x) - df[Solver_.Names.pcnv] = self.cc.value(x) - df[Solver_.Names.prad] = self.rc.value(x) - df[Solver_.Names.ppre] = self.pc.value(x) + df = self._steady_return_opt(return_err, return_power, x, z, err, df) return df @@ -539,24 +556,16 @@ def morgan(i: floatArray, tg: floatArray) -> floatArray: # format output df = pd.DataFrame({Solver_.Names.transit: x}) - - if return_err: - df["err"] = err - if return_temp or return_power: ts, tc = newtheader(x, y) ta = self.average(ts, tc) - if return_temp: df[Solver_.Names.tsurf] = ts df[Solver_.Names.tavg] = ta df[Solver_.Names.tcore] = tc - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(ta) - df[Solver_.Names.psol] = self.sh.value(ts) - df[Solver_.Names.pcnv] = self.cc.value(ts) - df[Solver_.Names.prad] = self.rc.value(ts) - df[Solver_.Names.ppre] = self.pc.value(ts) + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) return df From fb8342c80a808d6b71a894708ebe8458ae38bdb2 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:16:59 +0200 Subject: [PATCH 10/53] Reordered methods; added some type hints. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 136 ++++++++++---------- src/thermohl/solver/slv3t.py | 235 ++++++++++++++++++----------------- 2 files changed, 192 insertions(+), 179 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index a695a89..1eb8af4 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -96,6 +96,75 @@ def steady_temperature( return df + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + Imin: float = DP.imin, + Imax: float = DP.imax, + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_power: bool = True, + ) -> pd.DataFrame: + """Compute steady-state max intensity. + + Compute the maximum intensity that can be run in a conductor without + exceeding the temperature given in argument. + + Parameters + ---------- + T : float or numpy.ndarray + Maximum temperature. + Imin : float, optional + Lower bound for intensity. The default is 0. + Imax : float, optional + Upper bound for intensity. The default is 9999. + tol : float, optional + Tolerance for temperature error. The default is 1.0E-06. + maxiter : int, optional + Max number of iteration. The default is 64. + return_err : bool, optional + Return final error on intensity to check convergence. The default is False. + return_power : bool, optional + Return power term values. The default is True. + + Returns + ------- + pandas.DataFrame + A dataframe with maximum intensity and other results (depending on inputs) + in the columns. + + """ + + # save transit in arg + transit = self.args.I + + # solve with bisection + shape = self._min_shape() + T_ = T * np.ones(shape) + jh = ( + self.cc.value(T_) + + self.rc.value(T_) + + self.pc.value(T_) + - self.sh.value(T_) + ) + + def fun(i: floatArray) -> floatArrayLike: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + return self.jh.value(T_) - jh + + A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) + + # restore previous transit + self.args.I = transit + + # format output + df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) + df = self._steady_return_opt(return_err, return_power, T_, err, df) + + return df + def transient_temperature( self, time: floatArray = np.array([]), @@ -200,71 +269,4 @@ def transient_temperature( return dr - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - Imin: float = DP.imin, - Imax: float = DP.imax, - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_power: bool = True, - ) -> pd.DataFrame: - """Compute steady-state max intensity. - - Compute the maximum intensity that can be run in a conductor without - exceeding the temperature given in argument. - - Parameters - ---------- - T : float or numpy.ndarray - Maximum temperature. - Imin : float, optional - Lower bound for intensity. The default is 0. - Imax : float, optional - Upper bound for intensity. The default is 9999. - tol : float, optional - Tolerance for temperature error. The default is 1.0E-06. - maxiter : int, optional - Max number of iteration. The default is 64. - return_err : bool, optional - Return final error on intensity to check convergence. The default is False. - return_power : bool, optional - Return power term values. The default is True. - - Returns - ------- - pandas.DataFrame - A dataframe with maximum intensity and other results (depending on inputs) - in the columns. - - """ - - # save transit in arg - transit = self.args.I - - # solve with bisection - shape = self._min_shape() - T_ = T * np.ones(shape) - jh = ( - self.cc.value(T_) - + self.rc.value(T_) - + self.pc.value(T_) - - self.sh.value(T_) - ) - - def fun(i: floatArray) -> floatArrayLike: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - return self.jh.value(T_) - jh - - A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) - # restore previous transit - self.args.I = transit - - # format output - df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - df = self._steady_return_opt(return_err, return_power, T_, err, df) - - return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d8bd60a..40b29a1 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -209,6 +209,52 @@ def _steady_return_opt( return df + def _morgan_transient(self): + """Morgan coefficients for transient temperature.""" + c, D, d, ix = self.mgc + c1 = c / (2.0 * np.pi * self.args.l) + c2 = 0.5 * np.ones_like(c1) + a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) + c2[ix] = a / b + return c1, c2 + + def _transient_temperature_results( + self, + time: np.ndarray, + ts: np.ndarray, + ta: np.ndarray, + tc: np.ndarray, + return_power: bool, + n: int, + ): + """Format transient temperature results.""" + dr = { + Solver_.Names.time: time, + Solver_.Names.tsurf: ts, + Solver_.Names.tavg: ta, + Solver_.Names.tcore: tc, + } + + if return_power: + for power in Solver_.Names.powers(): + dr[power] = np.zeros_like(ts) + + for i in range(len(time)): + dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) + dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) + dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) + dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) + dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + + if n == 1: + keys = list(dr.keys()) + keys.remove(Solver_.Names.time) + for k in keys: + dr[k] = dr[k][:, 0] + + return dr + + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -270,41 +316,88 @@ def steady_temperature( return df - def _morgan_transient(self): - """Morgan coefficients for transient temperature.""" - c, D, d, ix = self.mgc - c1 = c / (2.0 * np.pi * self.args.l) - c2 = 0.5 * np.ones_like(c1) - a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) - c2[ix] = a / b - return c1, c2 + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + target: strListLike = "auto", + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_temp: bool = True, + return_power: bool = True, + ) -> pd.DataFrame: + """ + Compute the steady-state intensity for a given temperature profile. + Parameters: + ----------- + T : float or numpy.ndarray, optional + Initial temperature profile. Default is an empty numpy array. + target : str or List[str], optional + Target specification for the solver. Default is "auto". + tol : float, optional + Tolerance for the solver. Default is DP.tol. + maxiter : int, optional + Maximum number of iterations for the solver. Default is DP.maxiter. + return_err : bool, optional + If True, return the error in the output DataFrame. Default is False. + return_temp : bool, optional + If True, return the temperature profiles in the output DataFrame. Default is True. + return_power : bool, optional + If True, return the power profiles in the output DataFrame. Default is True. + Returns: + -------- + pd.DataFrame + DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. + """ - def _transient_temperature_results(self, time, ts, ta, tc, return_power, n): - dr = { - Solver_.Names.time: time, - Solver_.Names.tsurf: ts, - Solver_.Names.tavg: ta, - Solver_.Names.tcore: tc, - } + Tmax, newtheader = self._steady_intensity_header(T, target) - if return_power: - for power in Solver_.Names.powers(): - dr[power] = np.zeros_like(ts) + def balance(i: floatArray, tg: floatArray) -> floatArrayLike: + ts, tc = newtheader(i, tg) + return self.balance(ts, tc) - for i in range(len(time)): - dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) - dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) - dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) - dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) - dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + def morgan(i: floatArray, tg: floatArray) -> floatArray: + ts, tc = newtheader(i, tg) + return self.morgan(ts, tc) - if n == 1: - keys = list(dr.keys()) - keys.remove(Solver_.Names.time) - for k in keys: - dr[k] = dr[k][:, 0] + # solve system + s = Solver1T( + self.args.__dict__, + type(self.jh), + type(self.sh), + type(self.cc), + type(self.rc), + type(self.pc), + ) + r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) + x, y, cnt, err = quasi_newton_2d( + balance, + morgan, + r[Solver_.Names.transit].values, + Tmax, + relative_tolerance=tol, + max_iterations=maxiter, + delta_x=1.0e-03, + delta_y=1.0e-03, + ) + if np.max(err) > tol or cnt == maxiter: + print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - return dr + # format output + df = pd.DataFrame({Solver_.Names.transit: x}) + if return_temp or return_power: + ts, tc = newtheader(x, y) + ta = self.average(ts, tc) + if return_temp: + df[Solver_.Names.tsurf] = ts + df[Solver_.Names.tavg] = ta + df[Solver_.Names.tcore] = tc + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) + + return df def transient_temperature( self, @@ -487,85 +580,3 @@ def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: return Tmax, newtheader - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - target: strListLike = "auto", - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_temp: bool = True, - return_power: bool = True, - ) -> pd.DataFrame: - """ - Compute the steady-state intensity for a given temperature profile. - Parameters: - ----------- - T : float or numpy.ndarray, optional - Initial temperature profile. Default is an empty numpy array. - target : str or List[str], optional - Target specification for the solver. Default is "auto". - tol : float, optional - Tolerance for the solver. Default is DP.tol. - maxiter : int, optional - Maximum number of iterations for the solver. Default is DP.maxiter. - return_err : bool, optional - If True, return the error in the output DataFrame. Default is False. - return_temp : bool, optional - If True, return the temperature profiles in the output DataFrame. Default is True. - return_power : bool, optional - If True, return the power profiles in the output DataFrame. Default is True. - Returns: - -------- - pd.DataFrame - DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. - """ - - Tmax, newtheader = self._steady_intensity_header(T, target) - - def balance(i: floatArray, tg: floatArray) -> floatArrayLike: - ts, tc = newtheader(i, tg) - return self.balance(ts, tc) - - def morgan(i: floatArray, tg: floatArray) -> floatArray: - ts, tc = newtheader(i, tg) - return self.morgan(ts, tc) - - # solve system - s = Solver1T( - self.args.__dict__, - type(self.jh), - type(self.sh), - type(self.cc), - type(self.rc), - type(self.pc), - ) - r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) - x, y, cnt, err = quasi_newton_2d( - balance, - morgan, - r[Solver_.Names.transit].values, - Tmax, - relative_tolerance=tol, - max_iterations=maxiter, - delta_x=1.0e-03, - delta_y=1.0e-03, - ) - if np.max(err) > tol or cnt == maxiter: - print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - - # format output - df = pd.DataFrame({Solver_.Names.transit: x}) - if return_temp or return_power: - ts, tc = newtheader(x, y) - ta = self.average(ts, tc) - if return_temp: - df[Solver_.Names.tsurf] = ts - df[Solver_.Names.tavg] = ta - df[Solver_.Names.tcore] = tc - else: - ts = None - ta = None - df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) - - return df From 799f66eb70c6634629d86613959746e8f1035796 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:26:05 +0200 Subject: [PATCH 11/53] Reordered methods; moved check_target out of class scope. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 2 +- src/thermohl/solver/slv3t.py | 192 ++++++++++++++-------------- src/thermohl/solver/slv3t_legacy.py | 10 +- 3 files changed, 106 insertions(+), 98 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 1eb8af4..82278bb 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -12,7 +12,7 @@ import pandas as pd from thermohl import floatArrayLike, floatArray -from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.base import Solver as Solver_, Args from thermohl.solver.base import _DEFPARAM as DP from thermohl.solver.base import _set_dates, reshape from thermohl.utils import bisect_v diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 40b29a1..22fbc28 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -46,6 +46,61 @@ def _profile_bim_avg( return tc - (a / b) * (tc - ts) +def _check_target(target, d, max_len): + """ + Validates and processes the target temperature input for ampacity computations. + + Parameters: + target (str or list): The target temperature(s) to be validated. It can be: + - "auto": which sets the target to None. + - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + max_len (int): The expected length of the target list if target is a list. + + Returns: + numpy.ndarray: An array of target temperatures if the input is valid. + + Raises: + ValueError: If the target is a string not in the allowed values, or if the + target list length does not match max_len, or if any element in the target + list is not in the allowed values. + """ + + if target == "auto": + d_ = d * np.ones(max_len) + target_ = np.array( + [ + Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg + for i in range(max_len) + ] + ) + elif isinstance(target, str): + if target not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError( + f"Target temperature should be in " + f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" + f" got {target} instead." + ) + else: + target_ = np.array([target for _ in range(max_len)]) + else: + if len(target) != max_len: + raise ValueError() + for t in target: + if t not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError() + target_ = np.array(target) + return target_ + + class Solver3T(Solver_): def __init__( @@ -183,6 +238,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + # ========================================================================== + def _steady_return_opt( self, return_err: bool, @@ -209,6 +266,45 @@ def _steady_return_opt( return df + def _steady_intensity_header( + self, T: floatArrayLike, target: strListLike + ) -> Tuple[np.ndarray, Callable]: + """Format input for ampacity solver.""" + + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = _check_target(target, self.args.d, shape[0]) + + # pre-compute indexes + c, D, d, ix = self.mgc + a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) + + js = np.nonzero(target_ == Solver_.Names.surf)[0] + ja = np.nonzero(target_ == Solver_.Names.avg)[0] + jc = np.nonzero(target_ == Solver_.Names.core)[0] + jx = np.intersect1d(ix, ja) + + # get correct input for quasi-newton solver + def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + ts = np.ones_like(tg) * np.nan + tc = np.ones_like(tg) * np.nan + + ts[js] = Tmax[js] + tc[js] = tg[js] + + ts[ja] = tg[ja] + tc[ja] = 2 * Tmax[ja] - ts[ja] + tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) + + tc[jc] = Tmax[jc] + ts[jc] = tg[jc] + + return ts, tc + + return Tmax, newtheader + def _morgan_transient(self): """Morgan coefficients for transient temperature.""" c, D, d, ix = self.mgc @@ -254,6 +350,7 @@ def _transient_temperature_results( return dr + # ========================================================================== def steady_temperature( self, @@ -487,96 +584,3 @@ def transient_temperature( ts[i, :] = tc[i, :] - mrg return self._transient_temperature_results(time, ts, ta, tc, return_power, n) - - @staticmethod - def _check_target(target, d, max_len): - """ - Validates and processes the target temperature input. - - Parameters: - target (str or list): The target temperature(s) to be validated. It can be: - - "auto": which sets the target to None. - - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - max_len (int): The expected length of the target list if target is a list. - - Returns: - numpy.ndarray: An array of target temperatures if the input is valid. - - Raises: - ValueError: If the target is a string not in the allowed values, or if the target list length does not match max_len, or if any element in the target list is not in the allowed values. - """ - # check target - if target == "auto": - d_ = d * np.ones(max_len) - target_ = np.array( - [ - Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg - for i in range(max_len) - ] - ) - elif isinstance(target, str): - if target not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError( - f"Target temperature should be in " - f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" - f" got {target} instead." - ) - else: - target_ = np.array([target for _ in range(max_len)]) - else: - if len(target) != max_len: - raise ValueError() - for t in target: - if t not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError() - target_ = np.array(target) - return target_ - - def _steady_intensity_header( - self, T: floatArrayLike, target: strListLike - ) -> Tuple[np.ndarray, Callable]: - """Format input for ampacity solver.""" - - shape = self._min_shape() - Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) - - # pre-compute indexes - c, D, d, ix = self.mgc - a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) - - js = np.nonzero(target_ == Solver_.Names.surf)[0] - ja = np.nonzero(target_ == Solver_.Names.avg)[0] - jc = np.nonzero(target_ == Solver_.Names.core)[0] - jx = np.intersect1d(ix, ja) - - # get correct input for quasi-newton solver - def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - ts = np.ones_like(tg) * np.nan - tc = np.ones_like(tg) * np.nan - - ts[js] = Tmax[js] - tc[js] = tg[js] - - ts[ja] = tg[ja] - tc[ja] = 2 * Tmax[ja] - ts[ja] - tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) - - tc[jc] = Tmax[jc] - ts[jc] = tg[jc] - - return ts, tc - - return Tmax, newtheader - diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 561dfcb..53597f0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -11,8 +11,8 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, Args -from thermohl.solver.slv3t import Solver3T +from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.slv3t import Solver3T, _check_target class Solver3TL(Solver3T): @@ -86,6 +86,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: morgan_coefficient = self.mgc[0] return (tc - ts) - morgan_coefficient * self.joule(ts, tc) + # ========================================================================== + def _steady_intensity_header( self, T: floatArrayLike, target: strListLike ) -> Tuple[np.ndarray, Callable]: @@ -93,7 +95,7 @@ def _steady_intensity_header( shape = self._min_shape() Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) + target_ = _check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -125,6 +127,8 @@ def _morgan_transient(self): c2 = 0.5 * np.ones_like(c1) return c1, c2 + # ========================================================================== + def transient_temperature_legacy( self, time: floatArray = np.array([]), From 047fa2bb0a5ac175bb1de042b7de6b7bc351f359 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:29:59 +0200 Subject: [PATCH 12/53] Added dynamic parameter to transient temperature methods; factorized some code. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 31 ++++++++++++++++++ src/thermohl/solver/slv1t.py | 62 ++++++++++++------------------------ src/thermohl/solver/slv3t.py | 53 ++++++++++++------------------ 3 files changed, 71 insertions(+), 75 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 7034bfe..7eb4d39 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -260,6 +260,37 @@ def _min_shape(self) -> Tuple[int, ...]: shape = (1,) return shape + def _transient_process_dynamic( + self, time: np.ndarray, n: int, dynamic: dict = None + ): + """Code factorization for transient temperature computations. + + This methods prepare a dict with dynamic values to use in the + compute time loop. + """ + if len(time) < 2: + raise ValueError("The length of the time array must be at least 2.") + + # get month, day and hours in range with time + month, day, hour = _set_dates( + self.args.month, self.args.day, self.args.hour, time, n + ) + + # put dynamic values in a separate dict which will be used + # through the time loop + dynamic_ = { + "month": month, + "day": day, + "hour": hour, + } + + if dynamic is None: + dynamic = {} + for k, v in dynamic.items(): + dynamic_[k] = reshape(v, len(time), n) + + return dynamic_ + def update(self) -> None: self.args.extend() self.jh.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 82278bb..ca2d33c 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -169,6 +169,7 @@ def transient_temperature( self, time: floatArray = np.array([]), T0: Optional[float] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -194,49 +195,25 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError("The length of the time array must be at least 2.") - # get initial temperature - if T0 is None: - T0 = ( - self.args.Ta - if isinstance(self.args.Ta, numbers.Number) - else self.args.Ta[0] - ) - - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) - - # Two dicts, one (dc) with static quantities (with all elements of size n), the other (de) - # with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # save args + args = self.args.__dict__.copy() + + # initial conditions T = np.zeros((N, n)) + if T0 is None: + T0 = self.args.Ta T[0, :] = T0 - # main time loop - for i in range(1, len(time)): - for k, v in de.items(): + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): self.args[k] = v[i, :] self.update() T[i, :] = ( @@ -244,15 +221,15 @@ def transient_temperature( ) # save results - dr = dict(time=time, T=T) + dr = {Solver_.Names.time: time, Solver_.Names.temp: T} - # manage return dict 2 : powers + # add power to return dict if needed if return_power: for c in Solver_.Names.powers(): dr[c] = np.zeros_like(T) for i in range(N): - for k in de.keys(): - self.args[k] = de[k][i, :] + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() dr[Solver_.Names.pjle][i, :] = self.jh.value(T[i, :]) dr[Solver_.Names.psol][i, :] = self.sh.value(T[i, :]) @@ -260,13 +237,14 @@ def transient_temperature( dr[Solver_.Names.prad][i, :] = self.rc.value(T[i, :]) dr[Solver_.Names.ppre][i, :] = self.pc.value(T[i, :]) - # squeeze return values if n is 1 + # squeeze values in return dict (if n is 1) if n == 1: keys = list(dr.keys()) keys.remove(Solver_.Names.time) for k in keys: dr[k] = dr[k][:, 0] - return dr - + # restore args + self.args = Args(args) + return dr diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 22fbc28..023a547 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -501,6 +501,7 @@ def transient_temperature( time: floatArray = np.array([]), Ts0: Optional[floatArrayLike] = None, Tc0: Optional[floatArrayLike] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -531,51 +532,31 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError() - - # get initial temperature - Ts0 = Ts0 if Ts0 is not None else self.args.Ta - Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) - # Two dicts, one (dc) with static quantities (with all elements of size - # n), the other (de) with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # save args + args = self.args.__dict__.copy() # shortcuts for time-loop c1, c2 = self._morgan_transient() imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) - # main time loop - for i in range(1, len(time)): - for k in de.keys(): - self.args[k] = de[k][i, :] + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() bal = self.balance(ts[i - 1, :], tc[i - 1, :]) ta[i, :] = ta[i - 1, :] + (time[i] - time[i - 1]) * bal * imc @@ -583,4 +564,10 @@ def transient_temperature( tc[i, :] = ta[i, :] + c2 * mrg ts[i, :] = tc[i, :] - mrg - return self._transient_temperature_results(time, ts, ta, tc, return_power, n) + # get results + dr = self._transient_temperature_results(time, ts, ta, tc, return_power, n) + + # restore args + self.args = Args(args) + + return dr From bcc171cb22446fb92354a8862c08e6d55f5b31cd Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:30:25 +0200 Subject: [PATCH 13/53] Updated tests. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv3t_legacy.py | 6 +- test/functional_test/functional_test.py | 110 +++++++++++++++++++++++- test/unit/solver/test_slv1t.py | 14 +-- 3 files changed, 118 insertions(+), 12 deletions(-) diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 53597f0..b4171a0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -178,17 +178,19 @@ def transient_temperature_legacy( # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) dT = np.zeros((N, n)) - + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) dT[0, :] = tc[0, :] - ts[0, :] + # time loop for i in range(1, len(time)): bal = self.balance(ts[i - 1, :], tc[i - 1, :]) dti = time[i] - time[i - 1] diff --git a/test/functional_test/functional_test.py b/test/functional_test/functional_test.py index 19e927b..b9c164a 100644 --- a/test/functional_test/functional_test.py +++ b/test/functional_test/functional_test.py @@ -10,6 +10,7 @@ import numpy as np import pandas as pd +import thermohl.utils import yaml from typing import List, Dict @@ -67,13 +68,37 @@ def test_steady_temperature(): for d in scenario("temperature", "steady"): for _, e in d.items(): s = cner(scn2dict(e), heateq="3tl") - r = s.steady_temperature() + r = s.steady_temperature(return_power=False) assert np.allclose(r["t_surf"], e["T_surf"], atol=0.05) assert np.allclose(r["t_avg"], e["T_mean"], atol=0.05) assert np.allclose(r["t_core"], e["T_heart"], atol=0.05) +def test_steady_temperature_batch(): + din = [] + dref = [] + for d in scenario("temperature", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("T_surf", "T_mean", "T_heart") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + dref.rename( + columns={"T_surf": "t_surf", "T_mean": "t_avg", "T_heart": "t_core"}, + inplace=True, + ) + + s = cner(din, heateq="3tl") + r = s.steady_temperature(return_power=False) + assert np.allclose(r.values, dref.values, atol=0.05) + + def test_steady_ampacity(): for d in scenario("ampacity", "steady"): for _, e in d.items(): @@ -83,9 +108,29 @@ def test_steady_ampacity(): assert np.allclose(r["I"], e["I_max"], atol=0.05) +def test_steady_ampacity_batch(): + din = [] + dref = [] + for d in scenario("ampacity", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("Tmax_cable", "I_max") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + + s = cner(din, heateq="3tl") + r = s.steady_intensity(T=dref["Tmax_cable"], return_power=False) + assert np.allclose(r["I"].values, dref["I_max"].values, atol=0.05) + + def test_transient_temperature(): - atol = 0.5 + atol = 0.1 # this is hard-coded, maybe it should be put in the yaml file ... tau = 600.0 @@ -109,7 +154,7 @@ def test_transient_temperature(): rf = s.steady_temperature(Tsg=e["T_mean_final"], Tcg=e["T_mean_final"]) # time - time = np.arange(0.0, 1800.0, dt) + time = np.arange(0.0, 1800.0 + dt, dt) # transient temperature (linearized) rl = s.transient_temperature_legacy( @@ -128,3 +173,62 @@ def test_transient_temperature(): expected_temp = np.array(list(e[k1].values())) estimated_temp = np.interp(expected_time, rl["time"], rl[k2]) assert np.allclose(expected_temp, estimated_temp, atol=atol) + + +def test_transient_temperature_batch(): + + atol = 0.1 + + # this is hard-coded, maybe it should be put in the yaml file ... + tau = 600.0 + dt = 10.0 + minute = 60 + + din = [] + I0 = [] + dref = {"time": [], "t_surf": [], "t_avg": [], "t_core": [], "t_avg_final": []} + for d in scenario("temperature", "transient"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + I0.append(e["I0_cable"]) + for k1, k2 in zip( + ["T_surf_transient", "T_mean_transient", "T_heart_transient"], + ["t_surf", "t_avg", "t_core"], + ): + dref[k2].append(np.array(list(e[k1].values()))) + dref["t_avg_final"].append(e["T_mean_final"]) + din = pd.concat(din).reset_index(drop=True) + I0 = np.array(I0) + dref["time"] = np.array(list(e["T_surf_transient"].keys())) * minute * 1.0 + for k in ["t_surf", "t_avg", "t_core"]: + dref[k] = np.array(dref[k]).T + dref["t_avg_final"] = np.array(dref["t_avg_final"]) + + # solver + s = cner(din, heateq="3tl") + + # initial steady state + s.args["I"] = I0 + s.update() + ri = s.steady_temperature(return_power=False) + + # final steady state + s.args["I"] = din["I"] + s.update() + rf = s.steady_temperature(Tsg=dref["t_avg_final"], Tcg=dref["t_avg_final"]) + + # check final temp + assert np.allclose(dref["t_avg_final"], rf["t_avg"].values, atol=atol) + + # time + time = np.arange(0.0, 1800.0 + dt, dt) + + # transient temperature (linearized) + rl = s.transient_temperature_legacy( + time=time, Ts0=ri["t_surf"], Tc0=ri["t_core"], tau=tau + ) + + # filter expected times + ix = np.in1d(rl["time"], dref["time"]).nonzero()[0] + for k in ["t_surf", "t_avg", "t_core"]: + assert np.allclose(dref[k], rl[k][ix, :], atol=atol) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 4540faf..a413c43 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -88,9 +88,9 @@ def test_transient_temperature_default(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) def test_transient_temperature_with_initial_temp(solver): @@ -101,10 +101,10 @@ def test_transient_temperature_with_initial_temp(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) - assert result["T"][0] == T0 + assert len(result["t"]) == len(time) + assert result["t"][0] == T0 def test_transient_temperature_with_error(solver): @@ -114,9 +114,9 @@ def test_transient_temperature_with_error(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) assert Solver1T.Names.pjle in result assert Solver1T.Names.psol in result assert Solver1T.Names.pcnv in result From cfa9b1bc2cb8b38e0ac45809bf00389f64db13b8 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:33:28 +0200 Subject: [PATCH 14/53] Updated tests. Signed-off-by: Emmanuel Cieren --- test/test_shape.py | 21 +++++++++++++-------- test/test_solver3t.py | 1 + 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/test/test_shape.py b/test/test_shape.py index 20be2f4..243803a 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -29,14 +29,15 @@ def _ampargs(s: solver.Solver, t: pd.DataFrame): return a -def _traargs(s: solver.Solver, ds: pd.DataFrame, t, I): +def _traargs(s: solver.Solver, ds: pd.DataFrame, t: np.ndarray, dyn: dict): if isinstance(s, solver.Solver1T): - a = dict(time=t, T0=ds[solver.Solver.Names.temp].values) + a = dict(time=t, T0=ds[solver.Solver.Names.temp].values, dynamic=dyn) elif isinstance(s, solver.Solver3T): a = dict( time=t, Ts0=ds[solver.Solver.Names.tsurf].values, Tc0=ds[solver.Solver.Names.tcore].values, + dynamic=dyn, ) else: raise NotImplementedError @@ -104,13 +105,11 @@ def test_steady_1d_mix(): assert len(i) == n -def test_transient_0(): +def test_transient_no_dyn(): for s in _solvers(): t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) - ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, None) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) @@ -131,10 +130,16 @@ def test_transient_1(): s.update() t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) + d = { + "I": 199 * np.ones_like(t), + "Ta": 33.0, + "Pa": np.array(1.013e05), + "ws": 2.0 * np.ones(n), + "wa": np.ones((len(t), n)), + } ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, d) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) diff --git a/test/test_solver3t.py b/test/test_solver3t.py index 5dd4b73..fe2bc9f 100644 --- a/test/test_solver3t.py +++ b/test/test_solver3t.py @@ -86,6 +86,7 @@ def test_consistency(): target=t, return_err=True, return_power=True, + return_temp=True, tol=1.0e-09, maxiter=64, ) From 1e7ef39dfbfeec0b46e97bc9021fa03913b8a19f Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:02:38 +0200 Subject: [PATCH 15/53] Changed default value of srad parameter in all solar heating power to nan instead of None; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/power/cigre/solar_heating.py | 9 +++++---- src/thermohl/power/cner/solar_heating.py | 7 ++++--- src/thermohl/power/ieee/solar_heating.py | 7 ++++--- src/thermohl/power/olla/solar_heating.py | 9 +++++---- test/test_power_cner_integration.py | 2 +- 5 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/thermohl/power/cigre/solar_heating.py b/src/thermohl/power/cigre/solar_heating.py index d4482a8..eeed150 100644 --- a/src/thermohl/power/cigre/solar_heating.py +++ b/src/thermohl/power/cigre/solar_heating.py @@ -48,7 +48,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -74,8 +74,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx. + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns @@ -85,7 +86,7 @@ def __init__( """ self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = SolarHeating._solar_radiation( np.deg2rad(lat), np.deg2rad(azm), al, month, day, hour ) diff --git a/src/thermohl/power/cner/solar_heating.py b/src/thermohl/power/cner/solar_heating.py index 74b6e89..0e8a053 100644 --- a/src/thermohl/power/cner/solar_heating.py +++ b/src/thermohl/power/cner/solar_heating.py @@ -21,7 +21,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Build with args. @@ -45,8 +45,9 @@ def __init__( external diameter. alpha : np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/ieee/solar_heating.py b/src/thermohl/power/ieee/solar_heating.py index 24f00f9..0145670 100644 --- a/src/thermohl/power/ieee/solar_heating.py +++ b/src/thermohl/power/ieee/solar_heating.py @@ -24,7 +24,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -52,8 +52,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/olla/solar_heating.py b/src/thermohl/power/olla/solar_heating.py index 1df8d69..99bc5a7 100644 --- a/src/thermohl/power/olla/solar_heating.py +++ b/src/thermohl/power/olla/solar_heating.py @@ -24,7 +24,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -38,7 +38,7 @@ def __init__( lat : float or np.ndarray Latitude. alt : float or np.ndarray - Lltitude. + Longitude. azm : float or np.ndarray Azimuth. month : int or np.ndarray @@ -52,8 +52,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index d41ca0f..58b7ea6 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -220,6 +220,6 @@ def test_solar_heating(): alpha = 0.9 * ones p = np.array([34.9, 21.9357, 13.95, 21.9357, 21.9357]) - s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha, srad=None) + s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha) assert np.allclose(p, s.value(ones), 0.1) From 9001b13fa9b7f83c22c6a3f4b238056fdebd5afd Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:05:38 +0200 Subject: [PATCH 16/53] Changed default value of srad parameter in all base solar heating power term to nan instead of None Signed-off-by: Emmanuel Cieren --- src/thermohl/power/solar_heating.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/thermohl/power/solar_heating.py b/src/thermohl/power/solar_heating.py index b7a507a..2234a39 100644 --- a/src/thermohl/power/solar_heating.py +++ b/src/thermohl/power/solar_heating.py @@ -65,11 +65,11 @@ def __init__( D: floatArrayLike, alpha: floatArrayLike, est: _SRad, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = est(np.deg2rad(lat), alt, np.deg2rad(azm), tb, month, day, hour) else: self.srad = np.maximum(srad, 0.0) From 566c4b64d88c39de0ebfdeb94b6c87a0f79b1c50 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:06:25 +0200 Subject: [PATCH 17/53] Changed some var names to more readable, private values Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/__init__.py | 40 ++++++++++++++++----------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index b987948..36d2f30 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -9,10 +9,10 @@ from typing import Dict, Any, Optional, Union, Type -from thermohl.power import cigre as cigrep -from thermohl.power import cner as cnerp -from thermohl.power import ieee as ieeep -from thermohl.power import olla as ollap +from thermohl.power import cigre as _cigre +from thermohl.power import cner as _cner +from thermohl.power import ieee as _ieee +from thermohl.power import olla as _olla from thermohl.solver.base import Args, Solver from thermohl.solver.slv1d import Solver1D @@ -45,34 +45,34 @@ def _factory( if model == "cigre": return solver( dic, - cigrep.JouleHeating, - cigrep.SolarHeating, - cigrep.ConvectiveCooling, - cigrep.RadiativeCooling, + _cigre.JouleHeating, + _cigre.SolarHeating, + _cigre.ConvectiveCooling, + _cigre.RadiativeCooling, ) elif model == "ieee": return solver( dic, - ieeep.JouleHeating, - ieeep.SolarHeating, - ieeep.ConvectiveCooling, - ieeep.RadiativeCooling, + _ieee.JouleHeating, + _ieee.SolarHeating, + _ieee.ConvectiveCooling, + _ieee.RadiativeCooling, ) elif model == "olla": return solver( dic, - ollap.JouleHeating, - ollap.SolarHeating, - ollap.ConvectiveCooling, - ollap.RadiativeCooling, + _olla.JouleHeating, + _olla.SolarHeating, + _olla.ConvectiveCooling, + _olla.RadiativeCooling, ) elif model == "cner": return solver( dic, - cnerp.JouleHeating, - cnerp.SolarHeating, - cnerp.ConvectiveCooling, - cnerp.RadiativeCooling, + _cner.JouleHeating, + _cner.SolarHeating, + _cner.ConvectiveCooling, + _cner.RadiativeCooling, ) else: raise ValueError() From 4aa8fe0fe1b6e91c83fbcb316f24073c7176c9db Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:09:14 +0200 Subject: [PATCH 18/53] Added solar irradiance (srad) to Solver.args; updated comments Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index bc13416..1424a1f 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -62,9 +62,9 @@ def _set_default_values(self) -> None: self.alt = 0.0 # altitude (m) self.azm = 0.0 # azimuth (deg) - self.month = 3 # month number (1=Jan, 2=Feb, ...) + self.month = 3 # month number (1=Jan, 2=Feb ...) self.day = 21 # day of the month - self.hour = 12 # hour of the day (in [0, 23] range) + self.hour = 12.0 # hour of the day (in [0, 23] range) self.Ta = 15.0 # ambient temperature (C) self.Pa = 1.0e05 # ambient pressure (Pa) @@ -73,8 +73,8 @@ def _set_default_values(self) -> None: self.ws = 0.0 # wind speed (m.s**-1) self.wa = 90.0 # wind angle (deg, regarding north) self.al = 0.8 # albedo (1) - # coefficient for air pollution from 0 (clean) to 1 (polluted) - self.tb = 0.1 + self.tb = 0.1 # turbidity (coefficient for air pollution from 0 -clean- to 1 -polluted-) + self.srad = float("nan") # solar irradiance (in W.m**-2) self.I = 100.0 # transit intensity (A) From 07577e161876eeb1de64a8bb61c279d1f6810598 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:11:18 +0200 Subject: [PATCH 19/53] Renamed function with more explicit name; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/utils.py | 2 +- test/test_power_cner_integration.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index af82538..5fa21ac 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -101,7 +101,7 @@ def add_default_uncertainties(dat: dict, warning: bool = False) -> dict: return _dict_completion(dat, fil, check=False, warning=warning) -def df2dct(df: pd.DataFrame) -> dict: +def df2dict(df: pd.DataFrame) -> dict: """Convert a pandas.DataFrame to a dictionary. Would be an equivalent to df.to_dict(orient='numpy.ndarray') if it existed. diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index 58b7ea6..3f924cc 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -184,13 +184,13 @@ def test_compare_power(): ds = pd.concat(len(T) * (ds,)).reset_index(drop=True) T = np.concatenate([n * (t,) for t in T]) - from thermohl.utils import df2dct + from thermohl.utils import df2dict - d1 = df2dct(ds) + d1 = df2dict(ds) ds["wa"] = np.rad2deg( np.arcsin(np.sin(np.deg2rad(np.abs(ds["azm"] - ds["wa"]) % 180.0))) ) - d2 = df2dct(ds) + d2 = df2dict(ds) del (ds, n) pj = cner.JouleHeating(**d1) From fe8ed371a4371aec00f6f4b327d6b97ca323477c Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:53:12 +0200 Subject: [PATCH 20/53] Changed max_len/exten/compress methods in Solver class; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 98 +++++++++++++++-------------- src/thermohl/solver/slv1t.py | 11 +++- src/thermohl/solver/slv3t.py | 14 ++--- src/thermohl/solver/slv3t_legacy.py | 10 +-- test/test_shape.py | 2 +- test/unit/solver/test_base.py | 56 ++++++++--------- 6 files changed, 101 insertions(+), 90 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 1424a1f..4ec09f5 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -9,7 +9,7 @@ import datetime from abc import ABC, abstractmethod -from typing import Tuple, Type, Any, Optional, KeysView, Dict +from typing import Tuple, Type, Any, Optional, Dict, KeysView import numpy as np import pandas as pd @@ -53,6 +53,8 @@ def __init__(self, dic: Optional[dict[str, Any]] = None): for k in dic: if k in keys and dic[k] is not None: self[k] = dic[k] + # check for shape incompatibilities + _ = self.shape() def _set_default_values(self) -> None: """Set default values.""" @@ -115,53 +117,52 @@ def __getitem__(self, key: str) -> Any: def __setitem__(self, key: str, value: Any) -> None: self.__dict__[key] = value - def max_len(self) -> int: + def shape(self) -> Tuple[int, ...]: """ - Calculate the maximum length of the values in the dictionary. + Compute the maximum effective shape of values in current instance. - This method iterates over all keys in the dictionary and determines the maximum length - of the values associated with those keys. - If a value is not of a type that has a length, it is ignored. + Members of Args can be float of arrays. If arrays, they must be + one-diemnsional. Float and 1d array can coexist, but all arrays should + have the same shape/size. - Returns: - int: The maximum length of the values in the dictionary. If the dictionary is empty - or all values are of types that do not have a length, the method returns 1. - """ - n = 1 - for k in self.keys(): - try: - n = max(n, len(self[k])) - except TypeError: - pass - return n - - def extend_to_max_len(self) -> None: - """ - Extend all elements in the dictionary to the maximum length. + This method iterates over all keys in the instance's __dict__ and + computes the maximum length of the values associated with those keys. - This method iterates over all keys in the dictionary and checks if the - corresponding value is a numpy ndarray. If it is, it checks if its length - matches the maximum length obtained from the `max_len` method. - If the length matches, it creates a copy of the array. - If the length does not match or for non-ndarray values, it creates - a new numpy array of the maximum length, filled with the original value - and having the same data type. + If incompatible shapes are encountered, an exception is raised (ValueError). - Returns: - None """ - n = self.max_len() + shape_ = () for k in self.keys(): - if isinstance(self[k], np.ndarray): - t = self[k].dtype - c = len(self[k]) == n - else: - t = type(self[k]) - c = False - if c: - self[k] = self[k][:] + s = np.array(self[k]).shape + d = len(s) + er = f"Key {k} has a {s} shape when main shape is {shape_}" + if d == 0: + continue + if d == 1: + if shape_ == (): + shape_ = s + elif len(shape_) == 1: + if shape_ != s: + raise ValueError(er) + else: + raise ValueError(er) else: - self[k] = self[k] * np.ones((n,), dtype=t) + raise ValueError( + f"Key {k} has a {s} shape, only float and 1-dim arrays are accepted" + ) + return shape_ + + def extend(self, shape: Tuple[int, ...] = None) -> None: + # get shape + if shape is None: + shape = self.shape() + # complete if necessary + if len(shape) == 0: + shape = (1,) + # create a copy dict with scaled array + for k in self.keys(): + a = np.array(self[k]) + self[k] = a * np.ones(shape, dtype=a.dtype) def compress(self) -> None: """ @@ -172,10 +173,9 @@ def compress(self) -> None: None """ for k in self.keys(): - if isinstance(self[k], np.ndarray): - u = np.unique(self[k]) - if len(u) == 1: - self[k] = u[0] + u = np.unique(self[k]) + if len(u) == 1: + self[k] = u[0] class Solver(ABC): @@ -246,7 +246,7 @@ def __init__( """ self.args = Args(dic) - self.args.extend_to_max_len() + self.args.extend() self.jh = joule(**self.args.__dict__) self.sh = solar(**self.args.__dict__) self.cc = convective(**self.args.__dict__) @@ -254,8 +254,14 @@ def __init__( self.pc = precipitation(**self.args.__dict__) self.args.compress() + def _min_shape(self) -> Tuple[int, ...]: + shape = self.args.shape() + if shape == (): + shape = (1,) + return shape + def update(self) -> None: - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 03f9cfc..7fbdfbb 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -56,7 +56,12 @@ def steady_temperature( # solve with bisection T, err = bisect_v( - lambda x: -self.balance(x), Tmin, Tmax, (self.args.max_len(),), tol, maxiter + lambda x: -self.balance(x), + Tmin, + Tmax, + self._min_shape(), + tol=tol, + maxiter=maxiter, ) # format output @@ -101,7 +106,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError("The length of the time array must be at least 2.") @@ -222,7 +227,7 @@ def steady_intensity( transit = self.args.I # solve with bisection - shape = (self.args.max_len(),) + shape = self._min_shape() T_ = T * np.ones(shape) jh = ( self.cc.value(T_) diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index f038a0b..d04b790 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -78,7 +78,7 @@ def _morgan_coefficients( - i : numpy.ndarray[int] Indices where surface diameter `d_` is greater than 0. """ - c = 0.5 * np.ones((self.args.max_len(),)) + c = 0.5 * np.ones(self._min_shape()) D = self.args.D * np.ones_like(c) d = self.args.d * np.ones_like(c) i = np.nonzero(d > 0.0)[0] @@ -97,7 +97,7 @@ def update(self) -> None: Returns: None """ - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) @@ -215,7 +215,7 @@ def steady_temperature( """ # if no guess provided, use ambient temp - shape = (self.args.max_len(),) + shape = self._min_shape() Tsg = Tsg if Tsg is not None else 1.0 * self.args.Ta Tcg = Tcg if Tcg is not None else 1.5 * np.abs(self.args.Ta) Tsg_ = Tsg * np.ones(shape) @@ -322,7 +322,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() @@ -436,9 +436,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes c, D, d, ix = self.mgc diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 2da3a2f..561dfcb 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -46,7 +46,7 @@ def _morgan_coefficients(self) -> Tuple[floatArray, intArray]: UNIFORM_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 13 BIMETALLIC_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 21 - core_diameter_array = self.args.d * np.ones((self.args.max_len(),)) + core_diameter_array = self.args.d * np.ones(self._min_shape()) indices_non_zero_diameter = np.nonzero(core_diameter_array > 0.0)[0] heat_flux_coefficients = UNIFORM_CONDUCTOR_COEFFICIENT * np.ones_like( core_diameter_array @@ -91,9 +91,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -162,7 +162,7 @@ def transient_temperature_legacy( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() diff --git a/test/test_shape.py b/test/test_shape.py index 0c1f1ab..20be2f4 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -95,7 +95,7 @@ def test_steady_1d_mix(): n = 61 for s in _solvers(): s.args.Ta = np.linspace(-10, +50, n) - s.args.I = np.array([199.0]) + s.args.I = np.array(199.0) s.update() t = s.steady_temperature() a = _ampargs(s, t) diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 63eb841..9aba5d1 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -16,59 +16,59 @@ def test_max_len_with_mixed_types(): dic = { "lat": np.array([45.0, 46.0]), "lon": 10.0, - "alt": np.array([20.0, 25.0, 30.0]), + "alt": np.array([20.0, 25.0]), } args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 3 + assert result == (2,) -def test_max_len_with_ndarray(): +def test_shape_with_ndarray(): dic = {"lat": np.array([45.0, 46.0]), "lon": np.array([10.0, 11.0])} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 2 + assert result == (2,) -def test_max_len_with_scalar(): +def test_shape_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_empty_dict(): +def test_shape_with_empty_dict(): args = Args({}) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_varied_lengths(): +def test_shape_with_varied_lengths(): dic = { "lat": np.array([45.0, 46.0]), "lon": np.array([10.0]), "alt": np.array([20.0, 25.0, 30.0]), } - args = Args(dic) - - result = args.max_len() + try: + args = Args(dic) + assert False + except ValueError: + pass - assert result == 3 - -def test_extend_to_max_len_with_nd_array(): +def test_extend_with_nd_array(): dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -78,11 +78,11 @@ def test_extend_to_max_len_with_nd_array(): np.testing.assert_array_equal(args.lon, np.array([10.0, 10.0])) -def test_extend_to_max_len_with_scalar(): +def test_extend_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -92,11 +92,11 @@ def test_extend_to_max_len_with_scalar(): np.testing.assert_array_equal(args.lon, np.array([10.0])) -def test_extend_to_max_len_with_mixed_types(): - dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array([20.0])} +def test_extend_with_mixed_types(): + dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array(20.0)} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -109,10 +109,10 @@ def test_extend_to_max_len_with_mixed_types(): np.testing.assert_array_equal(args.alt, np.array([20.0, 20.0])) -def test_extend_to_max_len_with_empty_dict(): +def test_extend_with_empty_dict(): args = Args({}) - args.extend_to_max_len() + args.extend() for key in args.keys(): assert isinstance(args[key], (float, int, ndarray)) @@ -164,7 +164,7 @@ def test_compress_with_empty_dict(): args.compress() for key in args.keys(): - assert isinstance(args[key], (float, int, ndarray)) + assert isinstance(args[key], (float, int, np.integer, np.floating, ndarray)) if isinstance(args[key], ndarray): assert len(args[key]) == 1 From 1eabed4b5651578ceda21b5ee266621ef654b2dc Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:55:28 +0200 Subject: [PATCH 21/53] Fixed typo/error in set_times method; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 5 +---- test/unit/solver/test_base.py | 10 +++++----- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 4ec09f5..b1fe600 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -368,10 +368,7 @@ def _set_dates( td = np.array( [datetime.timedelta()] - + [ - datetime.timedelta(seconds=float(time[i] - time[i - 1])) - for i in range(1, N) - ] + + [datetime.timedelta(seconds=float(time[i] - time[0])) for i in range(1, N)] ) for j in range(n): diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 9aba5d1..089ed0c 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -256,7 +256,7 @@ def test_set_dates_single_day(): assert hours[1, 0] == 1.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 def test_set_dates_multiple_days(): @@ -279,7 +279,7 @@ def test_set_dates_multiple_days(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 2 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_months(): @@ -302,7 +302,7 @@ def test_set_dates_multiple_months(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_inputs(): @@ -328,7 +328,7 @@ def test_set_dates_multiple_inputs(): assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 assert months[0, 1] == 2 assert days[0, 1] == 2 @@ -340,4 +340,4 @@ def test_set_dates_multiple_inputs(): assert months[2, 1] == 2 assert days[2, 1] == 2 - assert hours[2, 1] == 13.0 + assert hours[2, 1] == 14.0 From eba5096d412e42363ff0a1df6f9614a95fc67cfa Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:58:08 +0200 Subject: [PATCH 22/53] Changed reshape function to be more restrictive; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 39 ++++++++++++++++++++--------------- test/unit/solver/test_base.py | 15 ++++++++++++-- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index b1fe600..7034bfe 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -291,12 +291,12 @@ def steady_intensity(self) -> pd.DataFrame: raise NotImplementedError -def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: +def reshape(input_var: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: """ Reshape the input array to the specified dimensions (nr, nc) if possible. Args: - input_array (numberArrayLike): Input array to be reshaped. + input (numberArrayLike): Variable to be reshaped. nb_row (int): Desired number of rows for the reshaped array. nb_columns (int): Desired number of columns for the reshaped array. @@ -305,24 +305,29 @@ def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numbe returns an array filled with the input_value repeated to fill the dimension (nb_row, nb_columns). Raises: - AttributeError: If the input_array has an invalid shape that cannot be reshaped. + ValueError: If the input has an invalid shape that cannot be reshaped. """ - reshaped_array = ndarray - try: - input_shape = input_array.shape - if len(input_shape) == 1: - if nb_row == input_shape[0]: - reshaped_array = np.column_stack(nb_columns * (input_array,)) - elif nb_columns == input_shape[0]: - reshaped_array = np.vstack(nb_row * (input_array,)) - elif len(input_shape) == 0: - raise AttributeError() - else: - reshaped_array = np.reshape(input_array, (nb_row, nb_columns)) - except AttributeError: + + input_array = np.array(input_var) + input_shape = input_array.shape + + msg = f"Input array has incompatible shape {input_shape} with specified number of rows ({nb_row}) and/or columns ({nb_columns})." + + if len(input_shape) == 0: reshaped_array = input_array * np.ones( - (nb_row, nb_columns), dtype=type(input_array) + (nb_row, nb_columns), dtype=input_array.dtype ) + elif len(input_shape) == 1: + if nb_row == input_shape[0]: + reshaped_array = np.column_stack(nb_columns * (input_array,)) + elif nb_columns == input_shape[0]: + reshaped_array = np.vstack(nb_row * (input_array,)) + else: + raise ValueError(msg) + elif input_shape == (nb_row, nb_columns): + return input_array + else: + raise ValueError(msg) return reshaped_array diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 089ed0c..f3822aa 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -225,8 +225,8 @@ def test_reshape_scalar_to_2d(): np.testing.assert_array_equal(result, expected) -def test_reshape_invalid_shape(): - array = np.array(0) # ([1.0, 2.0, 3.0]) +def test_reshape_another_scalar_to_2d(): + array = np.array(0) nb_row = 2 nb_columns = 2 expected = np.array([[0, 0], [0, 0]]) @@ -236,6 +236,17 @@ def test_reshape_invalid_shape(): np.testing.assert_array_equal(result, expected) +def test_reshape_invalid_shape(): + array = np.array([1.0, 2.0, 3.0]) + nb_row = 2 + nb_columns = 2 + try: + reshape(array, nb_row, nb_columns) + assert False + except ValueError: + pass + + def test_set_dates_single_day(): month = np.array([1]) day = np.array([1]) From ecb17beb717570fd1029eb573df0e021514665f9 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:07:46 +0200 Subject: [PATCH 23/53] Factorized some code for steady temperature methods. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 48 ++++++++++++++++++-------------- src/thermohl/solver/slv3t.py | 53 +++++++++++++++++++++--------------- 2 files changed, 59 insertions(+), 42 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 7fbdfbb..a695a89 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -19,6 +19,32 @@ class Solver1T(Solver_): + + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + T: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(T) + df[Solver_.Names.psol] = self.sh.value(T) + df[Solver_.Names.pcnv] = self.cc.value(T) + df[Solver_.Names.prad] = self.rc.value(T) + df[Solver_.Names.ppre] = self.pc.value(T) + + return df + def steady_temperature( self, Tmin: float = DP.tmin, @@ -66,16 +92,7 @@ def steady_temperature( # format output df = pd.DataFrame(data=T, columns=[Solver_.Names.temp]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T, err, df) return df @@ -248,15 +265,6 @@ def fun(i: floatArray) -> floatArrayLike: # format output df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T_, err, df) return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d04b790..d8bd60a 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -183,6 +183,32 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + Ts: np.ndarray, + Ta: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(Ta) + df[Solver_.Names.psol] = self.sh.value(Ts) + df[Solver_.Names.pcnv] = self.cc.value(Ts) + df[Solver_.Names.prad] = self.rc.value(Ts) + df[Solver_.Names.ppre] = self.pc.value(Ts) + + return df + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -240,16 +266,7 @@ def steady_temperature( df = pd.DataFrame( {Solver_.Names.tsurf: x, Solver_.Names.tavg: z, Solver_.Names.tcore: y} ) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.joule(x, y) - df[Solver_.Names.psol] = self.sh.value(x) - df[Solver_.Names.pcnv] = self.cc.value(x) - df[Solver_.Names.prad] = self.rc.value(x) - df[Solver_.Names.ppre] = self.pc.value(x) + df = self._steady_return_opt(return_err, return_power, x, z, err, df) return df @@ -539,24 +556,16 @@ def morgan(i: floatArray, tg: floatArray) -> floatArray: # format output df = pd.DataFrame({Solver_.Names.transit: x}) - - if return_err: - df["err"] = err - if return_temp or return_power: ts, tc = newtheader(x, y) ta = self.average(ts, tc) - if return_temp: df[Solver_.Names.tsurf] = ts df[Solver_.Names.tavg] = ta df[Solver_.Names.tcore] = tc - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(ta) - df[Solver_.Names.psol] = self.sh.value(ts) - df[Solver_.Names.pcnv] = self.cc.value(ts) - df[Solver_.Names.prad] = self.rc.value(ts) - df[Solver_.Names.ppre] = self.pc.value(ts) + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) return df From ff9c31e15d25349d17433df0dd5c4b9e6194acdd Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:16:59 +0200 Subject: [PATCH 24/53] Reordered methods; added some type hints. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 136 ++++++++++---------- src/thermohl/solver/slv3t.py | 235 ++++++++++++++++++----------------- 2 files changed, 192 insertions(+), 179 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index a695a89..1eb8af4 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -96,6 +96,75 @@ def steady_temperature( return df + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + Imin: float = DP.imin, + Imax: float = DP.imax, + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_power: bool = True, + ) -> pd.DataFrame: + """Compute steady-state max intensity. + + Compute the maximum intensity that can be run in a conductor without + exceeding the temperature given in argument. + + Parameters + ---------- + T : float or numpy.ndarray + Maximum temperature. + Imin : float, optional + Lower bound for intensity. The default is 0. + Imax : float, optional + Upper bound for intensity. The default is 9999. + tol : float, optional + Tolerance for temperature error. The default is 1.0E-06. + maxiter : int, optional + Max number of iteration. The default is 64. + return_err : bool, optional + Return final error on intensity to check convergence. The default is False. + return_power : bool, optional + Return power term values. The default is True. + + Returns + ------- + pandas.DataFrame + A dataframe with maximum intensity and other results (depending on inputs) + in the columns. + + """ + + # save transit in arg + transit = self.args.I + + # solve with bisection + shape = self._min_shape() + T_ = T * np.ones(shape) + jh = ( + self.cc.value(T_) + + self.rc.value(T_) + + self.pc.value(T_) + - self.sh.value(T_) + ) + + def fun(i: floatArray) -> floatArrayLike: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + return self.jh.value(T_) - jh + + A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) + + # restore previous transit + self.args.I = transit + + # format output + df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) + df = self._steady_return_opt(return_err, return_power, T_, err, df) + + return df + def transient_temperature( self, time: floatArray = np.array([]), @@ -200,71 +269,4 @@ def transient_temperature( return dr - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - Imin: float = DP.imin, - Imax: float = DP.imax, - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_power: bool = True, - ) -> pd.DataFrame: - """Compute steady-state max intensity. - - Compute the maximum intensity that can be run in a conductor without - exceeding the temperature given in argument. - - Parameters - ---------- - T : float or numpy.ndarray - Maximum temperature. - Imin : float, optional - Lower bound for intensity. The default is 0. - Imax : float, optional - Upper bound for intensity. The default is 9999. - tol : float, optional - Tolerance for temperature error. The default is 1.0E-06. - maxiter : int, optional - Max number of iteration. The default is 64. - return_err : bool, optional - Return final error on intensity to check convergence. The default is False. - return_power : bool, optional - Return power term values. The default is True. - - Returns - ------- - pandas.DataFrame - A dataframe with maximum intensity and other results (depending on inputs) - in the columns. - - """ - - # save transit in arg - transit = self.args.I - - # solve with bisection - shape = self._min_shape() - T_ = T * np.ones(shape) - jh = ( - self.cc.value(T_) - + self.rc.value(T_) - + self.pc.value(T_) - - self.sh.value(T_) - ) - - def fun(i: floatArray) -> floatArrayLike: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - return self.jh.value(T_) - jh - - A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) - # restore previous transit - self.args.I = transit - - # format output - df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - df = self._steady_return_opt(return_err, return_power, T_, err, df) - - return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d8bd60a..40b29a1 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -209,6 +209,52 @@ def _steady_return_opt( return df + def _morgan_transient(self): + """Morgan coefficients for transient temperature.""" + c, D, d, ix = self.mgc + c1 = c / (2.0 * np.pi * self.args.l) + c2 = 0.5 * np.ones_like(c1) + a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) + c2[ix] = a / b + return c1, c2 + + def _transient_temperature_results( + self, + time: np.ndarray, + ts: np.ndarray, + ta: np.ndarray, + tc: np.ndarray, + return_power: bool, + n: int, + ): + """Format transient temperature results.""" + dr = { + Solver_.Names.time: time, + Solver_.Names.tsurf: ts, + Solver_.Names.tavg: ta, + Solver_.Names.tcore: tc, + } + + if return_power: + for power in Solver_.Names.powers(): + dr[power] = np.zeros_like(ts) + + for i in range(len(time)): + dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) + dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) + dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) + dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) + dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + + if n == 1: + keys = list(dr.keys()) + keys.remove(Solver_.Names.time) + for k in keys: + dr[k] = dr[k][:, 0] + + return dr + + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -270,41 +316,88 @@ def steady_temperature( return df - def _morgan_transient(self): - """Morgan coefficients for transient temperature.""" - c, D, d, ix = self.mgc - c1 = c / (2.0 * np.pi * self.args.l) - c2 = 0.5 * np.ones_like(c1) - a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) - c2[ix] = a / b - return c1, c2 + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + target: strListLike = "auto", + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_temp: bool = True, + return_power: bool = True, + ) -> pd.DataFrame: + """ + Compute the steady-state intensity for a given temperature profile. + Parameters: + ----------- + T : float or numpy.ndarray, optional + Initial temperature profile. Default is an empty numpy array. + target : str or List[str], optional + Target specification for the solver. Default is "auto". + tol : float, optional + Tolerance for the solver. Default is DP.tol. + maxiter : int, optional + Maximum number of iterations for the solver. Default is DP.maxiter. + return_err : bool, optional + If True, return the error in the output DataFrame. Default is False. + return_temp : bool, optional + If True, return the temperature profiles in the output DataFrame. Default is True. + return_power : bool, optional + If True, return the power profiles in the output DataFrame. Default is True. + Returns: + -------- + pd.DataFrame + DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. + """ - def _transient_temperature_results(self, time, ts, ta, tc, return_power, n): - dr = { - Solver_.Names.time: time, - Solver_.Names.tsurf: ts, - Solver_.Names.tavg: ta, - Solver_.Names.tcore: tc, - } + Tmax, newtheader = self._steady_intensity_header(T, target) - if return_power: - for power in Solver_.Names.powers(): - dr[power] = np.zeros_like(ts) + def balance(i: floatArray, tg: floatArray) -> floatArrayLike: + ts, tc = newtheader(i, tg) + return self.balance(ts, tc) - for i in range(len(time)): - dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) - dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) - dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) - dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) - dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + def morgan(i: floatArray, tg: floatArray) -> floatArray: + ts, tc = newtheader(i, tg) + return self.morgan(ts, tc) - if n == 1: - keys = list(dr.keys()) - keys.remove(Solver_.Names.time) - for k in keys: - dr[k] = dr[k][:, 0] + # solve system + s = Solver1T( + self.args.__dict__, + type(self.jh), + type(self.sh), + type(self.cc), + type(self.rc), + type(self.pc), + ) + r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) + x, y, cnt, err = quasi_newton_2d( + balance, + morgan, + r[Solver_.Names.transit].values, + Tmax, + relative_tolerance=tol, + max_iterations=maxiter, + delta_x=1.0e-03, + delta_y=1.0e-03, + ) + if np.max(err) > tol or cnt == maxiter: + print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - return dr + # format output + df = pd.DataFrame({Solver_.Names.transit: x}) + if return_temp or return_power: + ts, tc = newtheader(x, y) + ta = self.average(ts, tc) + if return_temp: + df[Solver_.Names.tsurf] = ts + df[Solver_.Names.tavg] = ta + df[Solver_.Names.tcore] = tc + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) + + return df def transient_temperature( self, @@ -487,85 +580,3 @@ def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: return Tmax, newtheader - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - target: strListLike = "auto", - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_temp: bool = True, - return_power: bool = True, - ) -> pd.DataFrame: - """ - Compute the steady-state intensity for a given temperature profile. - Parameters: - ----------- - T : float or numpy.ndarray, optional - Initial temperature profile. Default is an empty numpy array. - target : str or List[str], optional - Target specification for the solver. Default is "auto". - tol : float, optional - Tolerance for the solver. Default is DP.tol. - maxiter : int, optional - Maximum number of iterations for the solver. Default is DP.maxiter. - return_err : bool, optional - If True, return the error in the output DataFrame. Default is False. - return_temp : bool, optional - If True, return the temperature profiles in the output DataFrame. Default is True. - return_power : bool, optional - If True, return the power profiles in the output DataFrame. Default is True. - Returns: - -------- - pd.DataFrame - DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. - """ - - Tmax, newtheader = self._steady_intensity_header(T, target) - - def balance(i: floatArray, tg: floatArray) -> floatArrayLike: - ts, tc = newtheader(i, tg) - return self.balance(ts, tc) - - def morgan(i: floatArray, tg: floatArray) -> floatArray: - ts, tc = newtheader(i, tg) - return self.morgan(ts, tc) - - # solve system - s = Solver1T( - self.args.__dict__, - type(self.jh), - type(self.sh), - type(self.cc), - type(self.rc), - type(self.pc), - ) - r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) - x, y, cnt, err = quasi_newton_2d( - balance, - morgan, - r[Solver_.Names.transit].values, - Tmax, - relative_tolerance=tol, - max_iterations=maxiter, - delta_x=1.0e-03, - delta_y=1.0e-03, - ) - if np.max(err) > tol or cnt == maxiter: - print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - - # format output - df = pd.DataFrame({Solver_.Names.transit: x}) - if return_temp or return_power: - ts, tc = newtheader(x, y) - ta = self.average(ts, tc) - if return_temp: - df[Solver_.Names.tsurf] = ts - df[Solver_.Names.tavg] = ta - df[Solver_.Names.tcore] = tc - else: - ts = None - ta = None - df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) - - return df From 5762a0f87baa047033ada096fabe883211041344 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:26:05 +0200 Subject: [PATCH 25/53] Reordered methods; moved check_target out of class scope. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 2 +- src/thermohl/solver/slv3t.py | 192 ++++++++++++++-------------- src/thermohl/solver/slv3t_legacy.py | 10 +- 3 files changed, 106 insertions(+), 98 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 1eb8af4..82278bb 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -12,7 +12,7 @@ import pandas as pd from thermohl import floatArrayLike, floatArray -from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.base import Solver as Solver_, Args from thermohl.solver.base import _DEFPARAM as DP from thermohl.solver.base import _set_dates, reshape from thermohl.utils import bisect_v diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 40b29a1..22fbc28 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -46,6 +46,61 @@ def _profile_bim_avg( return tc - (a / b) * (tc - ts) +def _check_target(target, d, max_len): + """ + Validates and processes the target temperature input for ampacity computations. + + Parameters: + target (str or list): The target temperature(s) to be validated. It can be: + - "auto": which sets the target to None. + - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + max_len (int): The expected length of the target list if target is a list. + + Returns: + numpy.ndarray: An array of target temperatures if the input is valid. + + Raises: + ValueError: If the target is a string not in the allowed values, or if the + target list length does not match max_len, or if any element in the target + list is not in the allowed values. + """ + + if target == "auto": + d_ = d * np.ones(max_len) + target_ = np.array( + [ + Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg + for i in range(max_len) + ] + ) + elif isinstance(target, str): + if target not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError( + f"Target temperature should be in " + f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" + f" got {target} instead." + ) + else: + target_ = np.array([target for _ in range(max_len)]) + else: + if len(target) != max_len: + raise ValueError() + for t in target: + if t not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError() + target_ = np.array(target) + return target_ + + class Solver3T(Solver_): def __init__( @@ -183,6 +238,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + # ========================================================================== + def _steady_return_opt( self, return_err: bool, @@ -209,6 +266,45 @@ def _steady_return_opt( return df + def _steady_intensity_header( + self, T: floatArrayLike, target: strListLike + ) -> Tuple[np.ndarray, Callable]: + """Format input for ampacity solver.""" + + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = _check_target(target, self.args.d, shape[0]) + + # pre-compute indexes + c, D, d, ix = self.mgc + a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) + + js = np.nonzero(target_ == Solver_.Names.surf)[0] + ja = np.nonzero(target_ == Solver_.Names.avg)[0] + jc = np.nonzero(target_ == Solver_.Names.core)[0] + jx = np.intersect1d(ix, ja) + + # get correct input for quasi-newton solver + def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + ts = np.ones_like(tg) * np.nan + tc = np.ones_like(tg) * np.nan + + ts[js] = Tmax[js] + tc[js] = tg[js] + + ts[ja] = tg[ja] + tc[ja] = 2 * Tmax[ja] - ts[ja] + tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) + + tc[jc] = Tmax[jc] + ts[jc] = tg[jc] + + return ts, tc + + return Tmax, newtheader + def _morgan_transient(self): """Morgan coefficients for transient temperature.""" c, D, d, ix = self.mgc @@ -254,6 +350,7 @@ def _transient_temperature_results( return dr + # ========================================================================== def steady_temperature( self, @@ -487,96 +584,3 @@ def transient_temperature( ts[i, :] = tc[i, :] - mrg return self._transient_temperature_results(time, ts, ta, tc, return_power, n) - - @staticmethod - def _check_target(target, d, max_len): - """ - Validates and processes the target temperature input. - - Parameters: - target (str or list): The target temperature(s) to be validated. It can be: - - "auto": which sets the target to None. - - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - max_len (int): The expected length of the target list if target is a list. - - Returns: - numpy.ndarray: An array of target temperatures if the input is valid. - - Raises: - ValueError: If the target is a string not in the allowed values, or if the target list length does not match max_len, or if any element in the target list is not in the allowed values. - """ - # check target - if target == "auto": - d_ = d * np.ones(max_len) - target_ = np.array( - [ - Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg - for i in range(max_len) - ] - ) - elif isinstance(target, str): - if target not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError( - f"Target temperature should be in " - f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" - f" got {target} instead." - ) - else: - target_ = np.array([target for _ in range(max_len)]) - else: - if len(target) != max_len: - raise ValueError() - for t in target: - if t not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError() - target_ = np.array(target) - return target_ - - def _steady_intensity_header( - self, T: floatArrayLike, target: strListLike - ) -> Tuple[np.ndarray, Callable]: - """Format input for ampacity solver.""" - - shape = self._min_shape() - Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) - - # pre-compute indexes - c, D, d, ix = self.mgc - a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) - - js = np.nonzero(target_ == Solver_.Names.surf)[0] - ja = np.nonzero(target_ == Solver_.Names.avg)[0] - jc = np.nonzero(target_ == Solver_.Names.core)[0] - jx = np.intersect1d(ix, ja) - - # get correct input for quasi-newton solver - def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - ts = np.ones_like(tg) * np.nan - tc = np.ones_like(tg) * np.nan - - ts[js] = Tmax[js] - tc[js] = tg[js] - - ts[ja] = tg[ja] - tc[ja] = 2 * Tmax[ja] - ts[ja] - tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) - - tc[jc] = Tmax[jc] - ts[jc] = tg[jc] - - return ts, tc - - return Tmax, newtheader - diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 561dfcb..53597f0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -11,8 +11,8 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, Args -from thermohl.solver.slv3t import Solver3T +from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.slv3t import Solver3T, _check_target class Solver3TL(Solver3T): @@ -86,6 +86,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: morgan_coefficient = self.mgc[0] return (tc - ts) - morgan_coefficient * self.joule(ts, tc) + # ========================================================================== + def _steady_intensity_header( self, T: floatArrayLike, target: strListLike ) -> Tuple[np.ndarray, Callable]: @@ -93,7 +95,7 @@ def _steady_intensity_header( shape = self._min_shape() Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) + target_ = _check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -125,6 +127,8 @@ def _morgan_transient(self): c2 = 0.5 * np.ones_like(c1) return c1, c2 + # ========================================================================== + def transient_temperature_legacy( self, time: floatArray = np.array([]), From f18f2f5d9aefaa98551c4533dc04f442db7e221f Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:29:59 +0200 Subject: [PATCH 26/53] Added dynamic parameter to transient temperature methods; factorized some code. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 31 ++++++++++++++++++ src/thermohl/solver/slv1t.py | 62 ++++++++++++------------------------ src/thermohl/solver/slv3t.py | 53 ++++++++++++------------------ 3 files changed, 71 insertions(+), 75 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 7034bfe..7eb4d39 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -260,6 +260,37 @@ def _min_shape(self) -> Tuple[int, ...]: shape = (1,) return shape + def _transient_process_dynamic( + self, time: np.ndarray, n: int, dynamic: dict = None + ): + """Code factorization for transient temperature computations. + + This methods prepare a dict with dynamic values to use in the + compute time loop. + """ + if len(time) < 2: + raise ValueError("The length of the time array must be at least 2.") + + # get month, day and hours in range with time + month, day, hour = _set_dates( + self.args.month, self.args.day, self.args.hour, time, n + ) + + # put dynamic values in a separate dict which will be used + # through the time loop + dynamic_ = { + "month": month, + "day": day, + "hour": hour, + } + + if dynamic is None: + dynamic = {} + for k, v in dynamic.items(): + dynamic_[k] = reshape(v, len(time), n) + + return dynamic_ + def update(self) -> None: self.args.extend() self.jh.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 82278bb..ca2d33c 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -169,6 +169,7 @@ def transient_temperature( self, time: floatArray = np.array([]), T0: Optional[float] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -194,49 +195,25 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError("The length of the time array must be at least 2.") - # get initial temperature - if T0 is None: - T0 = ( - self.args.Ta - if isinstance(self.args.Ta, numbers.Number) - else self.args.Ta[0] - ) - - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) - - # Two dicts, one (dc) with static quantities (with all elements of size n), the other (de) - # with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # save args + args = self.args.__dict__.copy() + + # initial conditions T = np.zeros((N, n)) + if T0 is None: + T0 = self.args.Ta T[0, :] = T0 - # main time loop - for i in range(1, len(time)): - for k, v in de.items(): + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): self.args[k] = v[i, :] self.update() T[i, :] = ( @@ -244,15 +221,15 @@ def transient_temperature( ) # save results - dr = dict(time=time, T=T) + dr = {Solver_.Names.time: time, Solver_.Names.temp: T} - # manage return dict 2 : powers + # add power to return dict if needed if return_power: for c in Solver_.Names.powers(): dr[c] = np.zeros_like(T) for i in range(N): - for k in de.keys(): - self.args[k] = de[k][i, :] + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() dr[Solver_.Names.pjle][i, :] = self.jh.value(T[i, :]) dr[Solver_.Names.psol][i, :] = self.sh.value(T[i, :]) @@ -260,13 +237,14 @@ def transient_temperature( dr[Solver_.Names.prad][i, :] = self.rc.value(T[i, :]) dr[Solver_.Names.ppre][i, :] = self.pc.value(T[i, :]) - # squeeze return values if n is 1 + # squeeze values in return dict (if n is 1) if n == 1: keys = list(dr.keys()) keys.remove(Solver_.Names.time) for k in keys: dr[k] = dr[k][:, 0] - return dr - + # restore args + self.args = Args(args) + return dr diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 22fbc28..023a547 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -501,6 +501,7 @@ def transient_temperature( time: floatArray = np.array([]), Ts0: Optional[floatArrayLike] = None, Tc0: Optional[floatArrayLike] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -531,51 +532,31 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError() - - # get initial temperature - Ts0 = Ts0 if Ts0 is not None else self.args.Ta - Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) - # Two dicts, one (dc) with static quantities (with all elements of size - # n), the other (de) with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # save args + args = self.args.__dict__.copy() # shortcuts for time-loop c1, c2 = self._morgan_transient() imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) - # main time loop - for i in range(1, len(time)): - for k in de.keys(): - self.args[k] = de[k][i, :] + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() bal = self.balance(ts[i - 1, :], tc[i - 1, :]) ta[i, :] = ta[i - 1, :] + (time[i] - time[i - 1]) * bal * imc @@ -583,4 +564,10 @@ def transient_temperature( tc[i, :] = ta[i, :] + c2 * mrg ts[i, :] = tc[i, :] - mrg - return self._transient_temperature_results(time, ts, ta, tc, return_power, n) + # get results + dr = self._transient_temperature_results(time, ts, ta, tc, return_power, n) + + # restore args + self.args = Args(args) + + return dr From 71da5b7e0d22ac8e8ad22a7c3972f2a107d805fe Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:30:25 +0200 Subject: [PATCH 27/53] Updated tests. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv3t_legacy.py | 6 +- test/functional_test/functional_test.py | 110 +++++++++++++++++++++++- test/unit/solver/test_slv1t.py | 14 +-- 3 files changed, 118 insertions(+), 12 deletions(-) diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 53597f0..b4171a0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -178,17 +178,19 @@ def transient_temperature_legacy( # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) dT = np.zeros((N, n)) - + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) dT[0, :] = tc[0, :] - ts[0, :] + # time loop for i in range(1, len(time)): bal = self.balance(ts[i - 1, :], tc[i - 1, :]) dti = time[i] - time[i - 1] diff --git a/test/functional_test/functional_test.py b/test/functional_test/functional_test.py index 19e927b..b9c164a 100644 --- a/test/functional_test/functional_test.py +++ b/test/functional_test/functional_test.py @@ -10,6 +10,7 @@ import numpy as np import pandas as pd +import thermohl.utils import yaml from typing import List, Dict @@ -67,13 +68,37 @@ def test_steady_temperature(): for d in scenario("temperature", "steady"): for _, e in d.items(): s = cner(scn2dict(e), heateq="3tl") - r = s.steady_temperature() + r = s.steady_temperature(return_power=False) assert np.allclose(r["t_surf"], e["T_surf"], atol=0.05) assert np.allclose(r["t_avg"], e["T_mean"], atol=0.05) assert np.allclose(r["t_core"], e["T_heart"], atol=0.05) +def test_steady_temperature_batch(): + din = [] + dref = [] + for d in scenario("temperature", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("T_surf", "T_mean", "T_heart") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + dref.rename( + columns={"T_surf": "t_surf", "T_mean": "t_avg", "T_heart": "t_core"}, + inplace=True, + ) + + s = cner(din, heateq="3tl") + r = s.steady_temperature(return_power=False) + assert np.allclose(r.values, dref.values, atol=0.05) + + def test_steady_ampacity(): for d in scenario("ampacity", "steady"): for _, e in d.items(): @@ -83,9 +108,29 @@ def test_steady_ampacity(): assert np.allclose(r["I"], e["I_max"], atol=0.05) +def test_steady_ampacity_batch(): + din = [] + dref = [] + for d in scenario("ampacity", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("Tmax_cable", "I_max") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + + s = cner(din, heateq="3tl") + r = s.steady_intensity(T=dref["Tmax_cable"], return_power=False) + assert np.allclose(r["I"].values, dref["I_max"].values, atol=0.05) + + def test_transient_temperature(): - atol = 0.5 + atol = 0.1 # this is hard-coded, maybe it should be put in the yaml file ... tau = 600.0 @@ -109,7 +154,7 @@ def test_transient_temperature(): rf = s.steady_temperature(Tsg=e["T_mean_final"], Tcg=e["T_mean_final"]) # time - time = np.arange(0.0, 1800.0, dt) + time = np.arange(0.0, 1800.0 + dt, dt) # transient temperature (linearized) rl = s.transient_temperature_legacy( @@ -128,3 +173,62 @@ def test_transient_temperature(): expected_temp = np.array(list(e[k1].values())) estimated_temp = np.interp(expected_time, rl["time"], rl[k2]) assert np.allclose(expected_temp, estimated_temp, atol=atol) + + +def test_transient_temperature_batch(): + + atol = 0.1 + + # this is hard-coded, maybe it should be put in the yaml file ... + tau = 600.0 + dt = 10.0 + minute = 60 + + din = [] + I0 = [] + dref = {"time": [], "t_surf": [], "t_avg": [], "t_core": [], "t_avg_final": []} + for d in scenario("temperature", "transient"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + I0.append(e["I0_cable"]) + for k1, k2 in zip( + ["T_surf_transient", "T_mean_transient", "T_heart_transient"], + ["t_surf", "t_avg", "t_core"], + ): + dref[k2].append(np.array(list(e[k1].values()))) + dref["t_avg_final"].append(e["T_mean_final"]) + din = pd.concat(din).reset_index(drop=True) + I0 = np.array(I0) + dref["time"] = np.array(list(e["T_surf_transient"].keys())) * minute * 1.0 + for k in ["t_surf", "t_avg", "t_core"]: + dref[k] = np.array(dref[k]).T + dref["t_avg_final"] = np.array(dref["t_avg_final"]) + + # solver + s = cner(din, heateq="3tl") + + # initial steady state + s.args["I"] = I0 + s.update() + ri = s.steady_temperature(return_power=False) + + # final steady state + s.args["I"] = din["I"] + s.update() + rf = s.steady_temperature(Tsg=dref["t_avg_final"], Tcg=dref["t_avg_final"]) + + # check final temp + assert np.allclose(dref["t_avg_final"], rf["t_avg"].values, atol=atol) + + # time + time = np.arange(0.0, 1800.0 + dt, dt) + + # transient temperature (linearized) + rl = s.transient_temperature_legacy( + time=time, Ts0=ri["t_surf"], Tc0=ri["t_core"], tau=tau + ) + + # filter expected times + ix = np.in1d(rl["time"], dref["time"]).nonzero()[0] + for k in ["t_surf", "t_avg", "t_core"]: + assert np.allclose(dref[k], rl[k][ix, :], atol=atol) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 3cbe419..f9626d9 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -87,9 +87,9 @@ def test_transient_temperature_default(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) def test_transient_temperature_with_initial_temp(solver): @@ -100,10 +100,10 @@ def test_transient_temperature_with_initial_temp(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) - assert result["T"][0] == T0 + assert len(result["t"]) == len(time) + assert result["t"][0] == T0 def test_transient_temperature_with_error(solver): @@ -113,9 +113,9 @@ def test_transient_temperature_with_error(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) assert Solver1T.Names.pjle in result assert Solver1T.Names.psol in result assert Solver1T.Names.pcnv in result From 59d6ef484acc467b91d37d8da913ee540c3ba280 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:33:28 +0200 Subject: [PATCH 28/53] Updated tests. Signed-off-by: Emmanuel Cieren --- test/test_shape.py | 21 +++++++++++++-------- test/test_solver3t.py | 1 + 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/test/test_shape.py b/test/test_shape.py index 20be2f4..243803a 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -29,14 +29,15 @@ def _ampargs(s: solver.Solver, t: pd.DataFrame): return a -def _traargs(s: solver.Solver, ds: pd.DataFrame, t, I): +def _traargs(s: solver.Solver, ds: pd.DataFrame, t: np.ndarray, dyn: dict): if isinstance(s, solver.Solver1T): - a = dict(time=t, T0=ds[solver.Solver.Names.temp].values) + a = dict(time=t, T0=ds[solver.Solver.Names.temp].values, dynamic=dyn) elif isinstance(s, solver.Solver3T): a = dict( time=t, Ts0=ds[solver.Solver.Names.tsurf].values, Tc0=ds[solver.Solver.Names.tcore].values, + dynamic=dyn, ) else: raise NotImplementedError @@ -104,13 +105,11 @@ def test_steady_1d_mix(): assert len(i) == n -def test_transient_0(): +def test_transient_no_dyn(): for s in _solvers(): t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) - ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, None) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) @@ -131,10 +130,16 @@ def test_transient_1(): s.update() t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) + d = { + "I": 199 * np.ones_like(t), + "Ta": 33.0, + "Pa": np.array(1.013e05), + "ws": 2.0 * np.ones(n), + "wa": np.ones((len(t), n)), + } ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, d) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) diff --git a/test/test_solver3t.py b/test/test_solver3t.py index 5dd4b73..fe2bc9f 100644 --- a/test/test_solver3t.py +++ b/test/test_solver3t.py @@ -86,6 +86,7 @@ def test_consistency(): target=t, return_err=True, return_power=True, + return_temp=True, tol=1.0e-09, maxiter=64, ) From c04de7d2b2996165239d377f1e7737ba762dc364 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 20:15:40 +0200 Subject: [PATCH 29/53] Updated sonar properties Signed-off-by: Emmanuel Cieren --- sonar-project.properties | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sonar-project.properties b/sonar-project.properties index 43ce9a2..c1ddfd3 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -6,8 +6,8 @@ # SPDX-License-Identifier: MPL-2.0 sonar.host.url=https://sonarcloud.io -sonar.organization=phlowers -sonar.projectKey=phlowers_thermohl +sonar.projectKey=ecieren_thermohl +sonar.organization=ecieren sonar.sources=src/ sonar.tests=test/ From 69922ab448fd631efa965bbae6e491906067bb0d Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:02:38 +0200 Subject: [PATCH 30/53] Changed default value of srad parameter in all solar heating power to nan instead of None; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/power/cigre/solar_heating.py | 9 +++++---- src/thermohl/power/cner/solar_heating.py | 7 ++++--- src/thermohl/power/ieee/solar_heating.py | 7 ++++--- src/thermohl/power/olla/solar_heating.py | 9 +++++---- test/test_power_cner_integration.py | 2 +- 5 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/thermohl/power/cigre/solar_heating.py b/src/thermohl/power/cigre/solar_heating.py index d4482a8..eeed150 100644 --- a/src/thermohl/power/cigre/solar_heating.py +++ b/src/thermohl/power/cigre/solar_heating.py @@ -48,7 +48,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -74,8 +74,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx. + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns @@ -85,7 +86,7 @@ def __init__( """ self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = SolarHeating._solar_radiation( np.deg2rad(lat), np.deg2rad(azm), al, month, day, hour ) diff --git a/src/thermohl/power/cner/solar_heating.py b/src/thermohl/power/cner/solar_heating.py index 74b6e89..0e8a053 100644 --- a/src/thermohl/power/cner/solar_heating.py +++ b/src/thermohl/power/cner/solar_heating.py @@ -21,7 +21,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Build with args. @@ -45,8 +45,9 @@ def __init__( external diameter. alpha : np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/ieee/solar_heating.py b/src/thermohl/power/ieee/solar_heating.py index 24f00f9..0145670 100644 --- a/src/thermohl/power/ieee/solar_heating.py +++ b/src/thermohl/power/ieee/solar_heating.py @@ -24,7 +24,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -52,8 +52,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/src/thermohl/power/olla/solar_heating.py b/src/thermohl/power/olla/solar_heating.py index 1df8d69..99bc5a7 100644 --- a/src/thermohl/power/olla/solar_heating.py +++ b/src/thermohl/power/olla/solar_heating.py @@ -24,7 +24,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Init with args. @@ -38,7 +38,7 @@ def __init__( lat : float or np.ndarray Latitude. alt : float or np.ndarray - Lltitude. + Longitude. azm : float or np.ndarray Azimuth. month : int or np.ndarray @@ -52,8 +52,9 @@ def __init__( external diameter. alpha : float or np.ndarray Solar absorption coefficient. - srad : xxx - xxx + srad : float or np.ndarray + Solar irradiance. Default is nan. If nan value is estimated using + all other input. Returns ------- diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index d41ca0f..58b7ea6 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -220,6 +220,6 @@ def test_solar_heating(): alpha = 0.9 * ones p = np.array([34.9, 21.9357, 13.95, 21.9357, 21.9357]) - s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha, srad=None) + s = cner.SolarHeating(lat, azm, month, day, hour, D, alpha) assert np.allclose(p, s.value(ones), 0.1) From 6e4ccf8115bd736ac0f7c2323d6aa56cab4ac9e6 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:05:38 +0200 Subject: [PATCH 31/53] Changed default value of srad parameter in all base solar heating power term to nan instead of None Signed-off-by: Emmanuel Cieren --- src/thermohl/power/solar_heating.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/thermohl/power/solar_heating.py b/src/thermohl/power/solar_heating.py index b7a507a..2234a39 100644 --- a/src/thermohl/power/solar_heating.py +++ b/src/thermohl/power/solar_heating.py @@ -65,11 +65,11 @@ def __init__( D: floatArrayLike, alpha: floatArrayLike, est: _SRad, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): self.alpha = alpha - if srad is None: + if np.all(np.isnan(srad)): self.srad = est(np.deg2rad(lat), alt, np.deg2rad(azm), tb, month, day, hour) else: self.srad = np.maximum(srad, 0.0) From 4f5ece325645074080cdf7a4d8e0e96cb98adc30 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:06:25 +0200 Subject: [PATCH 32/53] Changed some var names to more readable, private values Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/__init__.py | 40 ++++++++++++++++----------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index b987948..36d2f30 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -9,10 +9,10 @@ from typing import Dict, Any, Optional, Union, Type -from thermohl.power import cigre as cigrep -from thermohl.power import cner as cnerp -from thermohl.power import ieee as ieeep -from thermohl.power import olla as ollap +from thermohl.power import cigre as _cigre +from thermohl.power import cner as _cner +from thermohl.power import ieee as _ieee +from thermohl.power import olla as _olla from thermohl.solver.base import Args, Solver from thermohl.solver.slv1d import Solver1D @@ -45,34 +45,34 @@ def _factory( if model == "cigre": return solver( dic, - cigrep.JouleHeating, - cigrep.SolarHeating, - cigrep.ConvectiveCooling, - cigrep.RadiativeCooling, + _cigre.JouleHeating, + _cigre.SolarHeating, + _cigre.ConvectiveCooling, + _cigre.RadiativeCooling, ) elif model == "ieee": return solver( dic, - ieeep.JouleHeating, - ieeep.SolarHeating, - ieeep.ConvectiveCooling, - ieeep.RadiativeCooling, + _ieee.JouleHeating, + _ieee.SolarHeating, + _ieee.ConvectiveCooling, + _ieee.RadiativeCooling, ) elif model == "olla": return solver( dic, - ollap.JouleHeating, - ollap.SolarHeating, - ollap.ConvectiveCooling, - ollap.RadiativeCooling, + _olla.JouleHeating, + _olla.SolarHeating, + _olla.ConvectiveCooling, + _olla.RadiativeCooling, ) elif model == "cner": return solver( dic, - cnerp.JouleHeating, - cnerp.SolarHeating, - cnerp.ConvectiveCooling, - cnerp.RadiativeCooling, + _cner.JouleHeating, + _cner.SolarHeating, + _cner.ConvectiveCooling, + _cner.RadiativeCooling, ) else: raise ValueError() From 8966b3fa757153097f0ab1a3573e9d220f1e44fd Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:09:14 +0200 Subject: [PATCH 33/53] Added solar irradiance (srad) to Solver.args; updated comments Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index bc13416..1424a1f 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -62,9 +62,9 @@ def _set_default_values(self) -> None: self.alt = 0.0 # altitude (m) self.azm = 0.0 # azimuth (deg) - self.month = 3 # month number (1=Jan, 2=Feb, ...) + self.month = 3 # month number (1=Jan, 2=Feb ...) self.day = 21 # day of the month - self.hour = 12 # hour of the day (in [0, 23] range) + self.hour = 12.0 # hour of the day (in [0, 23] range) self.Ta = 15.0 # ambient temperature (C) self.Pa = 1.0e05 # ambient pressure (Pa) @@ -73,8 +73,8 @@ def _set_default_values(self) -> None: self.ws = 0.0 # wind speed (m.s**-1) self.wa = 90.0 # wind angle (deg, regarding north) self.al = 0.8 # albedo (1) - # coefficient for air pollution from 0 (clean) to 1 (polluted) - self.tb = 0.1 + self.tb = 0.1 # turbidity (coefficient for air pollution from 0 -clean- to 1 -polluted-) + self.srad = float("nan") # solar irradiance (in W.m**-2) self.I = 100.0 # transit intensity (A) From 5682553a9cab1265e2c35277559a392ac2c11c9a Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:11:18 +0200 Subject: [PATCH 34/53] Renamed function with more explicit name; updated test Signed-off-by: Emmanuel Cieren --- src/thermohl/utils.py | 2 +- test/test_power_cner_integration.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/thermohl/utils.py b/src/thermohl/utils.py index af82538..5fa21ac 100644 --- a/src/thermohl/utils.py +++ b/src/thermohl/utils.py @@ -101,7 +101,7 @@ def add_default_uncertainties(dat: dict, warning: bool = False) -> dict: return _dict_completion(dat, fil, check=False, warning=warning) -def df2dct(df: pd.DataFrame) -> dict: +def df2dict(df: pd.DataFrame) -> dict: """Convert a pandas.DataFrame to a dictionary. Would be an equivalent to df.to_dict(orient='numpy.ndarray') if it existed. diff --git a/test/test_power_cner_integration.py b/test/test_power_cner_integration.py index 58b7ea6..3f924cc 100644 --- a/test/test_power_cner_integration.py +++ b/test/test_power_cner_integration.py @@ -184,13 +184,13 @@ def test_compare_power(): ds = pd.concat(len(T) * (ds,)).reset_index(drop=True) T = np.concatenate([n * (t,) for t in T]) - from thermohl.utils import df2dct + from thermohl.utils import df2dict - d1 = df2dct(ds) + d1 = df2dict(ds) ds["wa"] = np.rad2deg( np.arcsin(np.sin(np.deg2rad(np.abs(ds["azm"] - ds["wa"]) % 180.0))) ) - d2 = df2dct(ds) + d2 = df2dict(ds) del (ds, n) pj = cner.JouleHeating(**d1) From 8c13dee35ad33073a3a883362a2a7f39e2a6d8fb Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:53:12 +0200 Subject: [PATCH 35/53] Changed max_len/exten/compress methods in Solver class; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 98 +++++++++++++++-------------- src/thermohl/solver/slv1t.py | 11 +++- src/thermohl/solver/slv3t.py | 14 ++--- src/thermohl/solver/slv3t_legacy.py | 10 +-- test/test_shape.py | 2 +- test/unit/solver/test_base.py | 56 ++++++++--------- 6 files changed, 101 insertions(+), 90 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 1424a1f..4ec09f5 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -9,7 +9,7 @@ import datetime from abc import ABC, abstractmethod -from typing import Tuple, Type, Any, Optional, KeysView, Dict +from typing import Tuple, Type, Any, Optional, Dict, KeysView import numpy as np import pandas as pd @@ -53,6 +53,8 @@ def __init__(self, dic: Optional[dict[str, Any]] = None): for k in dic: if k in keys and dic[k] is not None: self[k] = dic[k] + # check for shape incompatibilities + _ = self.shape() def _set_default_values(self) -> None: """Set default values.""" @@ -115,53 +117,52 @@ def __getitem__(self, key: str) -> Any: def __setitem__(self, key: str, value: Any) -> None: self.__dict__[key] = value - def max_len(self) -> int: + def shape(self) -> Tuple[int, ...]: """ - Calculate the maximum length of the values in the dictionary. + Compute the maximum effective shape of values in current instance. - This method iterates over all keys in the dictionary and determines the maximum length - of the values associated with those keys. - If a value is not of a type that has a length, it is ignored. + Members of Args can be float of arrays. If arrays, they must be + one-diemnsional. Float and 1d array can coexist, but all arrays should + have the same shape/size. - Returns: - int: The maximum length of the values in the dictionary. If the dictionary is empty - or all values are of types that do not have a length, the method returns 1. - """ - n = 1 - for k in self.keys(): - try: - n = max(n, len(self[k])) - except TypeError: - pass - return n - - def extend_to_max_len(self) -> None: - """ - Extend all elements in the dictionary to the maximum length. + This method iterates over all keys in the instance's __dict__ and + computes the maximum length of the values associated with those keys. - This method iterates over all keys in the dictionary and checks if the - corresponding value is a numpy ndarray. If it is, it checks if its length - matches the maximum length obtained from the `max_len` method. - If the length matches, it creates a copy of the array. - If the length does not match or for non-ndarray values, it creates - a new numpy array of the maximum length, filled with the original value - and having the same data type. + If incompatible shapes are encountered, an exception is raised (ValueError). - Returns: - None """ - n = self.max_len() + shape_ = () for k in self.keys(): - if isinstance(self[k], np.ndarray): - t = self[k].dtype - c = len(self[k]) == n - else: - t = type(self[k]) - c = False - if c: - self[k] = self[k][:] + s = np.array(self[k]).shape + d = len(s) + er = f"Key {k} has a {s} shape when main shape is {shape_}" + if d == 0: + continue + if d == 1: + if shape_ == (): + shape_ = s + elif len(shape_) == 1: + if shape_ != s: + raise ValueError(er) + else: + raise ValueError(er) else: - self[k] = self[k] * np.ones((n,), dtype=t) + raise ValueError( + f"Key {k} has a {s} shape, only float and 1-dim arrays are accepted" + ) + return shape_ + + def extend(self, shape: Tuple[int, ...] = None) -> None: + # get shape + if shape is None: + shape = self.shape() + # complete if necessary + if len(shape) == 0: + shape = (1,) + # create a copy dict with scaled array + for k in self.keys(): + a = np.array(self[k]) + self[k] = a * np.ones(shape, dtype=a.dtype) def compress(self) -> None: """ @@ -172,10 +173,9 @@ def compress(self) -> None: None """ for k in self.keys(): - if isinstance(self[k], np.ndarray): - u = np.unique(self[k]) - if len(u) == 1: - self[k] = u[0] + u = np.unique(self[k]) + if len(u) == 1: + self[k] = u[0] class Solver(ABC): @@ -246,7 +246,7 @@ def __init__( """ self.args = Args(dic) - self.args.extend_to_max_len() + self.args.extend() self.jh = joule(**self.args.__dict__) self.sh = solar(**self.args.__dict__) self.cc = convective(**self.args.__dict__) @@ -254,8 +254,14 @@ def __init__( self.pc = precipitation(**self.args.__dict__) self.args.compress() + def _min_shape(self) -> Tuple[int, ...]: + shape = self.args.shape() + if shape == (): + shape = (1,) + return shape + def update(self) -> None: - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 03f9cfc..7fbdfbb 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -56,7 +56,12 @@ def steady_temperature( # solve with bisection T, err = bisect_v( - lambda x: -self.balance(x), Tmin, Tmax, (self.args.max_len(),), tol, maxiter + lambda x: -self.balance(x), + Tmin, + Tmax, + self._min_shape(), + tol=tol, + maxiter=maxiter, ) # format output @@ -101,7 +106,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError("The length of the time array must be at least 2.") @@ -222,7 +227,7 @@ def steady_intensity( transit = self.args.I # solve with bisection - shape = (self.args.max_len(),) + shape = self._min_shape() T_ = T * np.ones(shape) jh = ( self.cc.value(T_) diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index f038a0b..d04b790 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -78,7 +78,7 @@ def _morgan_coefficients( - i : numpy.ndarray[int] Indices where surface diameter `d_` is greater than 0. """ - c = 0.5 * np.ones((self.args.max_len(),)) + c = 0.5 * np.ones(self._min_shape()) D = self.args.D * np.ones_like(c) d = self.args.d * np.ones_like(c) i = np.nonzero(d > 0.0)[0] @@ -97,7 +97,7 @@ def update(self) -> None: Returns: None """ - self.args.extend_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) @@ -215,7 +215,7 @@ def steady_temperature( """ # if no guess provided, use ambient temp - shape = (self.args.max_len(),) + shape = self._min_shape() Tsg = Tsg if Tsg is not None else 1.0 * self.args.Ta Tcg = Tcg if Tcg is not None else 1.5 * np.abs(self.args.Ta) Tsg_ = Tsg * np.ones(shape) @@ -322,7 +322,7 @@ def transient_temperature( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() @@ -436,9 +436,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes c, D, d, ix = self.mgc diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 2da3a2f..561dfcb 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -46,7 +46,7 @@ def _morgan_coefficients(self) -> Tuple[floatArray, intArray]: UNIFORM_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 13 BIMETALLIC_CONDUCTOR_COEFFICIENT: Final[float] = 1 / 21 - core_diameter_array = self.args.d * np.ones((self.args.max_len(),)) + core_diameter_array = self.args.d * np.ones(self._min_shape()) indices_non_zero_diameter = np.nonzero(core_diameter_array > 0.0)[0] heat_flux_coefficients = UNIFORM_CONDUCTOR_COEFFICIENT * np.ones_like( core_diameter_array @@ -91,9 +91,9 @@ def _steady_intensity_header( ) -> Tuple[np.ndarray, Callable]: """Format input for ampacity solver.""" - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = self._check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -162,7 +162,7 @@ def transient_temperature_legacy( """ # get sizes (n for input dict entries, N for time) - n = self.args.max_len() + n = self._min_shape()[0] N = len(time) if N < 2: raise ValueError() diff --git a/test/test_shape.py b/test/test_shape.py index 0c1f1ab..20be2f4 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -95,7 +95,7 @@ def test_steady_1d_mix(): n = 61 for s in _solvers(): s.args.Ta = np.linspace(-10, +50, n) - s.args.I = np.array([199.0]) + s.args.I = np.array(199.0) s.update() t = s.steady_temperature() a = _ampargs(s, t) diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 63eb841..9aba5d1 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -16,59 +16,59 @@ def test_max_len_with_mixed_types(): dic = { "lat": np.array([45.0, 46.0]), "lon": 10.0, - "alt": np.array([20.0, 25.0, 30.0]), + "alt": np.array([20.0, 25.0]), } args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 3 + assert result == (2,) -def test_max_len_with_ndarray(): +def test_shape_with_ndarray(): dic = {"lat": np.array([45.0, 46.0]), "lon": np.array([10.0, 11.0])} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 2 + assert result == (2,) -def test_max_len_with_scalar(): +def test_shape_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_empty_dict(): +def test_shape_with_empty_dict(): args = Args({}) - result = args.max_len() + result = args.shape() - assert result == 1 + assert result == () -def test_max_len_with_varied_lengths(): +def test_shape_with_varied_lengths(): dic = { "lat": np.array([45.0, 46.0]), "lon": np.array([10.0]), "alt": np.array([20.0, 25.0, 30.0]), } - args = Args(dic) - - result = args.max_len() + try: + args = Args(dic) + assert False + except ValueError: + pass - assert result == 3 - -def test_extend_to_max_len_with_nd_array(): +def test_extend_with_nd_array(): dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -78,11 +78,11 @@ def test_extend_to_max_len_with_nd_array(): np.testing.assert_array_equal(args.lon, np.array([10.0, 10.0])) -def test_extend_to_max_len_with_scalar(): +def test_extend_with_scalar(): dic = {"lat": 45.0, "lon": 10.0} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -92,11 +92,11 @@ def test_extend_to_max_len_with_scalar(): np.testing.assert_array_equal(args.lon, np.array([10.0])) -def test_extend_to_max_len_with_mixed_types(): - dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array([20.0])} +def test_extend_with_mixed_types(): + dic = {"lat": np.array([45.0, 46.0]), "lon": 10.0, "alt": np.array(20.0)} args = Args(dic) - args.extend_to_max_len() + args.extend() assert isinstance(args.lat, ndarray) assert isinstance(args.lon, ndarray) @@ -109,10 +109,10 @@ def test_extend_to_max_len_with_mixed_types(): np.testing.assert_array_equal(args.alt, np.array([20.0, 20.0])) -def test_extend_to_max_len_with_empty_dict(): +def test_extend_with_empty_dict(): args = Args({}) - args.extend_to_max_len() + args.extend() for key in args.keys(): assert isinstance(args[key], (float, int, ndarray)) @@ -164,7 +164,7 @@ def test_compress_with_empty_dict(): args.compress() for key in args.keys(): - assert isinstance(args[key], (float, int, ndarray)) + assert isinstance(args[key], (float, int, np.integer, np.floating, ndarray)) if isinstance(args[key], ndarray): assert len(args[key]) == 1 From a6a74c588c39432fa4df370daaa3acd8f79ab785 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:55:28 +0200 Subject: [PATCH 36/53] Fixed typo/error in set_times method; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 5 +---- test/unit/solver/test_base.py | 10 +++++----- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 4ec09f5..b1fe600 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -368,10 +368,7 @@ def _set_dates( td = np.array( [datetime.timedelta()] - + [ - datetime.timedelta(seconds=float(time[i] - time[i - 1])) - for i in range(1, N) - ] + + [datetime.timedelta(seconds=float(time[i] - time[0])) for i in range(1, N)] ) for j in range(n): diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 9aba5d1..089ed0c 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -256,7 +256,7 @@ def test_set_dates_single_day(): assert hours[1, 0] == 1.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 def test_set_dates_multiple_days(): @@ -279,7 +279,7 @@ def test_set_dates_multiple_days(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 2 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_months(): @@ -302,7 +302,7 @@ def test_set_dates_multiple_months(): assert hours[1, 0] == 0.0 assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 0.0 + assert hours[2, 0] == 1.0 def test_set_dates_multiple_inputs(): @@ -328,7 +328,7 @@ def test_set_dates_multiple_inputs(): assert months[2, 0] == 1 assert days[2, 0] == 1 - assert hours[2, 0] == 1.0 + assert hours[2, 0] == 2.0 assert months[0, 1] == 2 assert days[0, 1] == 2 @@ -340,4 +340,4 @@ def test_set_dates_multiple_inputs(): assert months[2, 1] == 2 assert days[2, 1] == 2 - assert hours[2, 1] == 13.0 + assert hours[2, 1] == 14.0 From 60c0ad73fc418d684dc96395d0a5e060f742287d Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:58:08 +0200 Subject: [PATCH 37/53] Changed reshape function to be more restrictive; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 39 ++++++++++++++++++++--------------- test/unit/solver/test_base.py | 15 ++++++++++++-- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index b1fe600..7034bfe 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -291,12 +291,12 @@ def steady_intensity(self) -> pd.DataFrame: raise NotImplementedError -def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: +def reshape(input_var: numberArrayLike, nb_row: int, nb_columns: int) -> numberArray: """ Reshape the input array to the specified dimensions (nr, nc) if possible. Args: - input_array (numberArrayLike): Input array to be reshaped. + input (numberArrayLike): Variable to be reshaped. nb_row (int): Desired number of rows for the reshaped array. nb_columns (int): Desired number of columns for the reshaped array. @@ -305,24 +305,29 @@ def reshape(input_array: numberArrayLike, nb_row: int, nb_columns: int) -> numbe returns an array filled with the input_value repeated to fill the dimension (nb_row, nb_columns). Raises: - AttributeError: If the input_array has an invalid shape that cannot be reshaped. + ValueError: If the input has an invalid shape that cannot be reshaped. """ - reshaped_array = ndarray - try: - input_shape = input_array.shape - if len(input_shape) == 1: - if nb_row == input_shape[0]: - reshaped_array = np.column_stack(nb_columns * (input_array,)) - elif nb_columns == input_shape[0]: - reshaped_array = np.vstack(nb_row * (input_array,)) - elif len(input_shape) == 0: - raise AttributeError() - else: - reshaped_array = np.reshape(input_array, (nb_row, nb_columns)) - except AttributeError: + + input_array = np.array(input_var) + input_shape = input_array.shape + + msg = f"Input array has incompatible shape {input_shape} with specified number of rows ({nb_row}) and/or columns ({nb_columns})." + + if len(input_shape) == 0: reshaped_array = input_array * np.ones( - (nb_row, nb_columns), dtype=type(input_array) + (nb_row, nb_columns), dtype=input_array.dtype ) + elif len(input_shape) == 1: + if nb_row == input_shape[0]: + reshaped_array = np.column_stack(nb_columns * (input_array,)) + elif nb_columns == input_shape[0]: + reshaped_array = np.vstack(nb_row * (input_array,)) + else: + raise ValueError(msg) + elif input_shape == (nb_row, nb_columns): + return input_array + else: + raise ValueError(msg) return reshaped_array diff --git a/test/unit/solver/test_base.py b/test/unit/solver/test_base.py index 089ed0c..f3822aa 100644 --- a/test/unit/solver/test_base.py +++ b/test/unit/solver/test_base.py @@ -225,8 +225,8 @@ def test_reshape_scalar_to_2d(): np.testing.assert_array_equal(result, expected) -def test_reshape_invalid_shape(): - array = np.array(0) # ([1.0, 2.0, 3.0]) +def test_reshape_another_scalar_to_2d(): + array = np.array(0) nb_row = 2 nb_columns = 2 expected = np.array([[0, 0], [0, 0]]) @@ -236,6 +236,17 @@ def test_reshape_invalid_shape(): np.testing.assert_array_equal(result, expected) +def test_reshape_invalid_shape(): + array = np.array([1.0, 2.0, 3.0]) + nb_row = 2 + nb_columns = 2 + try: + reshape(array, nb_row, nb_columns) + assert False + except ValueError: + pass + + def test_set_dates_single_day(): month = np.array([1]) day = np.array([1]) From bdc1d612d623bda86c288eaf32815325e6e37f88 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:07:46 +0200 Subject: [PATCH 38/53] Factorized some code for steady temperature methods. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 48 ++++++++++++++++++-------------- src/thermohl/solver/slv3t.py | 53 +++++++++++++++++++++--------------- 2 files changed, 59 insertions(+), 42 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 7fbdfbb..a695a89 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -19,6 +19,32 @@ class Solver1T(Solver_): + + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + T: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(T) + df[Solver_.Names.psol] = self.sh.value(T) + df[Solver_.Names.pcnv] = self.cc.value(T) + df[Solver_.Names.prad] = self.rc.value(T) + df[Solver_.Names.ppre] = self.pc.value(T) + + return df + def steady_temperature( self, Tmin: float = DP.tmin, @@ -66,16 +92,7 @@ def steady_temperature( # format output df = pd.DataFrame(data=T, columns=[Solver_.Names.temp]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T, err, df) return df @@ -248,15 +265,6 @@ def fun(i: floatArray) -> floatArrayLike: # format output df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(T) - df[Solver_.Names.psol] = self.sh.value(T) - df[Solver_.Names.pcnv] = self.cc.value(T) - df[Solver_.Names.prad] = self.rc.value(T) - df[Solver_.Names.ppre] = self.pc.value(T) + df = self._steady_return_opt(return_err, return_power, T_, err, df) return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d04b790..d8bd60a 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -183,6 +183,32 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + def _steady_return_opt( + self, + return_err: bool, + return_power: bool, + Ts: np.ndarray, + Ta: np.ndarray, + err: np.ndarray, + df: pd.DataFrame, + ): + """Add error and/or power values to pd.Dataframe returned in + steady_temperature and steady_intensity methods.""" + + # add convergence error if asked + if return_err: + df[Solver_.Names.err] = err + + # add power values if asked + if return_power: + df[Solver_.Names.pjle] = self.jh.value(Ta) + df[Solver_.Names.psol] = self.sh.value(Ts) + df[Solver_.Names.pcnv] = self.cc.value(Ts) + df[Solver_.Names.prad] = self.rc.value(Ts) + df[Solver_.Names.ppre] = self.pc.value(Ts) + + return df + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -240,16 +266,7 @@ def steady_temperature( df = pd.DataFrame( {Solver_.Names.tsurf: x, Solver_.Names.tavg: z, Solver_.Names.tcore: y} ) - - if return_err: - df[Solver_.Names.err] = err - - if return_power: - df[Solver_.Names.pjle] = self.joule(x, y) - df[Solver_.Names.psol] = self.sh.value(x) - df[Solver_.Names.pcnv] = self.cc.value(x) - df[Solver_.Names.prad] = self.rc.value(x) - df[Solver_.Names.ppre] = self.pc.value(x) + df = self._steady_return_opt(return_err, return_power, x, z, err, df) return df @@ -539,24 +556,16 @@ def morgan(i: floatArray, tg: floatArray) -> floatArray: # format output df = pd.DataFrame({Solver_.Names.transit: x}) - - if return_err: - df["err"] = err - if return_temp or return_power: ts, tc = newtheader(x, y) ta = self.average(ts, tc) - if return_temp: df[Solver_.Names.tsurf] = ts df[Solver_.Names.tavg] = ta df[Solver_.Names.tcore] = tc - - if return_power: - df[Solver_.Names.pjle] = self.jh.value(ta) - df[Solver_.Names.psol] = self.sh.value(ts) - df[Solver_.Names.pcnv] = self.cc.value(ts) - df[Solver_.Names.prad] = self.rc.value(ts) - df[Solver_.Names.ppre] = self.pc.value(ts) + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) return df From 56c8137eecb72c01ce9ece48e3b2b1a228bcc9ae Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:16:59 +0200 Subject: [PATCH 39/53] Reordered methods; added some type hints. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 136 ++++++++++---------- src/thermohl/solver/slv3t.py | 235 ++++++++++++++++++----------------- 2 files changed, 192 insertions(+), 179 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index a695a89..1eb8af4 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -96,6 +96,75 @@ def steady_temperature( return df + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + Imin: float = DP.imin, + Imax: float = DP.imax, + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_power: bool = True, + ) -> pd.DataFrame: + """Compute steady-state max intensity. + + Compute the maximum intensity that can be run in a conductor without + exceeding the temperature given in argument. + + Parameters + ---------- + T : float or numpy.ndarray + Maximum temperature. + Imin : float, optional + Lower bound for intensity. The default is 0. + Imax : float, optional + Upper bound for intensity. The default is 9999. + tol : float, optional + Tolerance for temperature error. The default is 1.0E-06. + maxiter : int, optional + Max number of iteration. The default is 64. + return_err : bool, optional + Return final error on intensity to check convergence. The default is False. + return_power : bool, optional + Return power term values. The default is True. + + Returns + ------- + pandas.DataFrame + A dataframe with maximum intensity and other results (depending on inputs) + in the columns. + + """ + + # save transit in arg + transit = self.args.I + + # solve with bisection + shape = self._min_shape() + T_ = T * np.ones(shape) + jh = ( + self.cc.value(T_) + + self.rc.value(T_) + + self.pc.value(T_) + - self.sh.value(T_) + ) + + def fun(i: floatArray) -> floatArrayLike: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + return self.jh.value(T_) - jh + + A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) + + # restore previous transit + self.args.I = transit + + # format output + df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) + df = self._steady_return_opt(return_err, return_power, T_, err, df) + + return df + def transient_temperature( self, time: floatArray = np.array([]), @@ -200,71 +269,4 @@ def transient_temperature( return dr - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - Imin: float = DP.imin, - Imax: float = DP.imax, - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_power: bool = True, - ) -> pd.DataFrame: - """Compute steady-state max intensity. - - Compute the maximum intensity that can be run in a conductor without - exceeding the temperature given in argument. - - Parameters - ---------- - T : float or numpy.ndarray - Maximum temperature. - Imin : float, optional - Lower bound for intensity. The default is 0. - Imax : float, optional - Upper bound for intensity. The default is 9999. - tol : float, optional - Tolerance for temperature error. The default is 1.0E-06. - maxiter : int, optional - Max number of iteration. The default is 64. - return_err : bool, optional - Return final error on intensity to check convergence. The default is False. - return_power : bool, optional - Return power term values. The default is True. - - Returns - ------- - pandas.DataFrame - A dataframe with maximum intensity and other results (depending on inputs) - in the columns. - - """ - - # save transit in arg - transit = self.args.I - - # solve with bisection - shape = self._min_shape() - T_ = T * np.ones(shape) - jh = ( - self.cc.value(T_) - + self.rc.value(T_) - + self.pc.value(T_) - - self.sh.value(T_) - ) - - def fun(i: floatArray) -> floatArrayLike: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - return self.jh.value(T_) - jh - - A, err = bisect_v(fun, Imin, Imax, shape, tol, maxiter) - # restore previous transit - self.args.I = transit - - # format output - df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) - df = self._steady_return_opt(return_err, return_power, T_, err, df) - - return df diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index d8bd60a..40b29a1 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -209,6 +209,52 @@ def _steady_return_opt( return df + def _morgan_transient(self): + """Morgan coefficients for transient temperature.""" + c, D, d, ix = self.mgc + c1 = c / (2.0 * np.pi * self.args.l) + c2 = 0.5 * np.ones_like(c1) + a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) + c2[ix] = a / b + return c1, c2 + + def _transient_temperature_results( + self, + time: np.ndarray, + ts: np.ndarray, + ta: np.ndarray, + tc: np.ndarray, + return_power: bool, + n: int, + ): + """Format transient temperature results.""" + dr = { + Solver_.Names.time: time, + Solver_.Names.tsurf: ts, + Solver_.Names.tavg: ta, + Solver_.Names.tcore: tc, + } + + if return_power: + for power in Solver_.Names.powers(): + dr[power] = np.zeros_like(ts) + + for i in range(len(time)): + dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) + dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) + dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) + dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) + dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + + if n == 1: + keys = list(dr.keys()) + keys.remove(Solver_.Names.time) + for k in keys: + dr[k] = dr[k][:, 0] + + return dr + + def steady_temperature( self, Tsg: Optional[floatArrayLike] = None, @@ -270,41 +316,88 @@ def steady_temperature( return df - def _morgan_transient(self): - """Morgan coefficients for transient temperature.""" - c, D, d, ix = self.mgc - c1 = c / (2.0 * np.pi * self.args.l) - c2 = 0.5 * np.ones_like(c1) - a, b = _profile_bim_avg_coeffs(0.5 * d[ix], 0.5 * D[ix]) - c2[ix] = a / b - return c1, c2 + def steady_intensity( + self, + T: floatArrayLike = np.array([]), + target: strListLike = "auto", + tol: float = DP.tol, + maxiter: int = DP.maxiter, + return_err: bool = False, + return_temp: bool = True, + return_power: bool = True, + ) -> pd.DataFrame: + """ + Compute the steady-state intensity for a given temperature profile. + Parameters: + ----------- + T : float or numpy.ndarray, optional + Initial temperature profile. Default is an empty numpy array. + target : str or List[str], optional + Target specification for the solver. Default is "auto". + tol : float, optional + Tolerance for the solver. Default is DP.tol. + maxiter : int, optional + Maximum number of iterations for the solver. Default is DP.maxiter. + return_err : bool, optional + If True, return the error in the output DataFrame. Default is False. + return_temp : bool, optional + If True, return the temperature profiles in the output DataFrame. Default is True. + return_power : bool, optional + If True, return the power profiles in the output DataFrame. Default is True. + Returns: + -------- + pd.DataFrame + DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. + """ - def _transient_temperature_results(self, time, ts, ta, tc, return_power, n): - dr = { - Solver_.Names.time: time, - Solver_.Names.tsurf: ts, - Solver_.Names.tavg: ta, - Solver_.Names.tcore: tc, - } + Tmax, newtheader = self._steady_intensity_header(T, target) - if return_power: - for power in Solver_.Names.powers(): - dr[power] = np.zeros_like(ts) + def balance(i: floatArray, tg: floatArray) -> floatArrayLike: + ts, tc = newtheader(i, tg) + return self.balance(ts, tc) - for i in range(len(time)): - dr[Solver_.Names.pjle][i, :] = self.joule(ts[i, :], tc[i, :]) - dr[Solver_.Names.psol][i, :] = self.sh.value(ts[i, :]) - dr[Solver_.Names.pcnv][i, :] = self.cc.value(ts[i, :]) - dr[Solver_.Names.prad][i, :] = self.rc.value(ts[i, :]) - dr[Solver_.Names.ppre][i, :] = self.pc.value(ts[i, :]) + def morgan(i: floatArray, tg: floatArray) -> floatArray: + ts, tc = newtheader(i, tg) + return self.morgan(ts, tc) - if n == 1: - keys = list(dr.keys()) - keys.remove(Solver_.Names.time) - for k in keys: - dr[k] = dr[k][:, 0] + # solve system + s = Solver1T( + self.args.__dict__, + type(self.jh), + type(self.sh), + type(self.cc), + type(self.rc), + type(self.pc), + ) + r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) + x, y, cnt, err = quasi_newton_2d( + balance, + morgan, + r[Solver_.Names.transit].values, + Tmax, + relative_tolerance=tol, + max_iterations=maxiter, + delta_x=1.0e-03, + delta_y=1.0e-03, + ) + if np.max(err) > tol or cnt == maxiter: + print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - return dr + # format output + df = pd.DataFrame({Solver_.Names.transit: x}) + if return_temp or return_power: + ts, tc = newtheader(x, y) + ta = self.average(ts, tc) + if return_temp: + df[Solver_.Names.tsurf] = ts + df[Solver_.Names.tavg] = ta + df[Solver_.Names.tcore] = tc + else: + ts = None + ta = None + df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) + + return df def transient_temperature( self, @@ -487,85 +580,3 @@ def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: return Tmax, newtheader - def steady_intensity( - self, - T: floatArrayLike = np.array([]), - target: strListLike = "auto", - tol: float = DP.tol, - maxiter: int = DP.maxiter, - return_err: bool = False, - return_temp: bool = True, - return_power: bool = True, - ) -> pd.DataFrame: - """ - Compute the steady-state intensity for a given temperature profile. - Parameters: - ----------- - T : float or numpy.ndarray, optional - Initial temperature profile. Default is an empty numpy array. - target : str or List[str], optional - Target specification for the solver. Default is "auto". - tol : float, optional - Tolerance for the solver. Default is DP.tol. - maxiter : int, optional - Maximum number of iterations for the solver. Default is DP.maxiter. - return_err : bool, optional - If True, return the error in the output DataFrame. Default is False. - return_temp : bool, optional - If True, return the temperature profiles in the output DataFrame. Default is True. - return_power : bool, optional - If True, return the power profiles in the output DataFrame. Default is True. - Returns: - -------- - pd.DataFrame - DataFrame containing the steady-state intensity and optionally the error, temperature profiles, and power profiles. - """ - - Tmax, newtheader = self._steady_intensity_header(T, target) - - def balance(i: floatArray, tg: floatArray) -> floatArrayLike: - ts, tc = newtheader(i, tg) - return self.balance(ts, tc) - - def morgan(i: floatArray, tg: floatArray) -> floatArray: - ts, tc = newtheader(i, tg) - return self.morgan(ts, tc) - - # solve system - s = Solver1T( - self.args.__dict__, - type(self.jh), - type(self.sh), - type(self.cc), - type(self.rc), - type(self.pc), - ) - r = s.steady_intensity(Tmax, tol=1.0, maxiter=8, return_power=False) - x, y, cnt, err = quasi_newton_2d( - balance, - morgan, - r[Solver_.Names.transit].values, - Tmax, - relative_tolerance=tol, - max_iterations=maxiter, - delta_x=1.0e-03, - delta_y=1.0e-03, - ) - if np.max(err) > tol or cnt == maxiter: - print(f"rstat_analytic max err is {np.max(err):.3E} in {cnt:d} iterations") - - # format output - df = pd.DataFrame({Solver_.Names.transit: x}) - if return_temp or return_power: - ts, tc = newtheader(x, y) - ta = self.average(ts, tc) - if return_temp: - df[Solver_.Names.tsurf] = ts - df[Solver_.Names.tavg] = ta - df[Solver_.Names.tcore] = tc - else: - ts = None - ta = None - df = self._steady_return_opt(return_err, return_power, ts, ta, err, df) - - return df From cc098ec15741e05921d6bf5e948396a0f1ccd74e Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:26:05 +0200 Subject: [PATCH 40/53] Reordered methods; moved check_target out of class scope. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv1t.py | 2 +- src/thermohl/solver/slv3t.py | 192 ++++++++++++++-------------- src/thermohl/solver/slv3t_legacy.py | 10 +- 3 files changed, 106 insertions(+), 98 deletions(-) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 1eb8af4..82278bb 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -12,7 +12,7 @@ import pandas as pd from thermohl import floatArrayLike, floatArray -from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.base import Solver as Solver_, Args from thermohl.solver.base import _DEFPARAM as DP from thermohl.solver.base import _set_dates, reshape from thermohl.utils import bisect_v diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 40b29a1..22fbc28 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -46,6 +46,61 @@ def _profile_bim_avg( return tc - (a / b) * (tc - ts) +def _check_target(target, d, max_len): + """ + Validates and processes the target temperature input for ampacity computations. + + Parameters: + target (str or list): The target temperature(s) to be validated. It can be: + - "auto": which sets the target to None. + - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. + max_len (int): The expected length of the target list if target is a list. + + Returns: + numpy.ndarray: An array of target temperatures if the input is valid. + + Raises: + ValueError: If the target is a string not in the allowed values, or if the + target list length does not match max_len, or if any element in the target + list is not in the allowed values. + """ + + if target == "auto": + d_ = d * np.ones(max_len) + target_ = np.array( + [ + Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg + for i in range(max_len) + ] + ) + elif isinstance(target, str): + if target not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError( + f"Target temperature should be in " + f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" + f" got {target} instead." + ) + else: + target_ = np.array([target for _ in range(max_len)]) + else: + if len(target) != max_len: + raise ValueError() + for t in target: + if t not in [ + Solver_.Names.surf, + Solver_.Names.avg, + Solver_.Names.core, + ]: + raise ValueError() + target_ = np.array(target) + return target_ + + class Solver3T(Solver_): def __init__( @@ -183,6 +238,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + # ========================================================================== + def _steady_return_opt( self, return_err: bool, @@ -209,6 +266,45 @@ def _steady_return_opt( return df + def _steady_intensity_header( + self, T: floatArrayLike, target: strListLike + ) -> Tuple[np.ndarray, Callable]: + """Format input for ampacity solver.""" + + shape = self._min_shape() + Tmax = T * np.ones(shape) + target_ = _check_target(target, self.args.d, shape[0]) + + # pre-compute indexes + c, D, d, ix = self.mgc + a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) + + js = np.nonzero(target_ == Solver_.Names.surf)[0] + ja = np.nonzero(target_ == Solver_.Names.avg)[0] + jc = np.nonzero(target_ == Solver_.Names.core)[0] + jx = np.intersect1d(ix, ja) + + # get correct input for quasi-newton solver + def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: + self.args.I = i + self.jh.__init__(**self.args.__dict__) + ts = np.ones_like(tg) * np.nan + tc = np.ones_like(tg) * np.nan + + ts[js] = Tmax[js] + tc[js] = tg[js] + + ts[ja] = tg[ja] + tc[ja] = 2 * Tmax[ja] - ts[ja] + tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) + + tc[jc] = Tmax[jc] + ts[jc] = tg[jc] + + return ts, tc + + return Tmax, newtheader + def _morgan_transient(self): """Morgan coefficients for transient temperature.""" c, D, d, ix = self.mgc @@ -254,6 +350,7 @@ def _transient_temperature_results( return dr + # ========================================================================== def steady_temperature( self, @@ -487,96 +584,3 @@ def transient_temperature( ts[i, :] = tc[i, :] - mrg return self._transient_temperature_results(time, ts, ta, tc, return_power, n) - - @staticmethod - def _check_target(target, d, max_len): - """ - Validates and processes the target temperature input. - - Parameters: - target (str or list): The target temperature(s) to be validated. It can be: - - "auto": which sets the target to None. - - A string: which should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - - A list of strings: where each string should be one of Solver_.Names.surf, Solver_.Names.avg, or Solver_.Names.core. - max_len (int): The expected length of the target list if target is a list. - - Returns: - numpy.ndarray: An array of target temperatures if the input is valid. - - Raises: - ValueError: If the target is a string not in the allowed values, or if the target list length does not match max_len, or if any element in the target list is not in the allowed values. - """ - # check target - if target == "auto": - d_ = d * np.ones(max_len) - target_ = np.array( - [ - Solver_.Names.core if d_[i] > 0.0 else Solver_.Names.avg - for i in range(max_len) - ] - ) - elif isinstance(target, str): - if target not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError( - f"Target temperature should be in " - f"{[Solver_.Names.surf, Solver_.Names.avg, Solver_.Names.core]};" - f" got {target} instead." - ) - else: - target_ = np.array([target for _ in range(max_len)]) - else: - if len(target) != max_len: - raise ValueError() - for t in target: - if t not in [ - Solver_.Names.surf, - Solver_.Names.avg, - Solver_.Names.core, - ]: - raise ValueError() - target_ = np.array(target) - return target_ - - def _steady_intensity_header( - self, T: floatArrayLike, target: strListLike - ) -> Tuple[np.ndarray, Callable]: - """Format input for ampacity solver.""" - - shape = self._min_shape() - Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) - - # pre-compute indexes - c, D, d, ix = self.mgc - a, b = _profile_bim_avg_coeffs(0.5 * d, 0.5 * D) - - js = np.nonzero(target_ == Solver_.Names.surf)[0] - ja = np.nonzero(target_ == Solver_.Names.avg)[0] - jc = np.nonzero(target_ == Solver_.Names.core)[0] - jx = np.intersect1d(ix, ja) - - # get correct input for quasi-newton solver - def newtheader(i: floatArray, tg: floatArray) -> Tuple[floatArray, floatArray]: - self.args.I = i - self.jh.__init__(**self.args.__dict__) - ts = np.ones_like(tg) * np.nan - tc = np.ones_like(tg) * np.nan - - ts[js] = Tmax[js] - tc[js] = tg[js] - - ts[ja] = tg[ja] - tc[ja] = 2 * Tmax[ja] - ts[ja] - tc[jx] = (b[jx] * Tmax[jx] - a[jx] * ts[jx]) / (b[jx] - a[jx]) - - tc[jc] = Tmax[jc] - ts[jc] = tg[jc] - - return ts, tc - - return Tmax, newtheader - diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 561dfcb..53597f0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -11,8 +11,8 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, Args -from thermohl.solver.slv3t import Solver3T +from thermohl.solver.base import Solver as Solver_ +from thermohl.solver.slv3t import Solver3T, _check_target class Solver3TL(Solver3T): @@ -86,6 +86,8 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: morgan_coefficient = self.mgc[0] return (tc - ts) - morgan_coefficient * self.joule(ts, tc) + # ========================================================================== + def _steady_intensity_header( self, T: floatArrayLike, target: strListLike ) -> Tuple[np.ndarray, Callable]: @@ -93,7 +95,7 @@ def _steady_intensity_header( shape = self._min_shape() Tmax = T * np.ones(shape) - target_ = self._check_target(target, self.args.d, shape[0]) + target_ = _check_target(target, self.args.d, shape[0]) # pre-compute indexes surface_indices = np.nonzero(target_ == Solver_.Names.surf)[0] @@ -125,6 +127,8 @@ def _morgan_transient(self): c2 = 0.5 * np.ones_like(c1) return c1, c2 + # ========================================================================== + def transient_temperature_legacy( self, time: floatArray = np.array([]), From 6d98323cf19975bd2409963742943e500d068ea2 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:29:59 +0200 Subject: [PATCH 41/53] Added dynamic parameter to transient temperature methods; factorized some code. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 31 ++++++++++++++++++ src/thermohl/solver/slv1t.py | 62 ++++++++++++------------------------ src/thermohl/solver/slv3t.py | 53 ++++++++++++------------------ 3 files changed, 71 insertions(+), 75 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 7034bfe..7eb4d39 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -260,6 +260,37 @@ def _min_shape(self) -> Tuple[int, ...]: shape = (1,) return shape + def _transient_process_dynamic( + self, time: np.ndarray, n: int, dynamic: dict = None + ): + """Code factorization for transient temperature computations. + + This methods prepare a dict with dynamic values to use in the + compute time loop. + """ + if len(time) < 2: + raise ValueError("The length of the time array must be at least 2.") + + # get month, day and hours in range with time + month, day, hour = _set_dates( + self.args.month, self.args.day, self.args.hour, time, n + ) + + # put dynamic values in a separate dict which will be used + # through the time loop + dynamic_ = { + "month": month, + "day": day, + "hour": hour, + } + + if dynamic is None: + dynamic = {} + for k, v in dynamic.items(): + dynamic_[k] = reshape(v, len(time), n) + + return dynamic_ + def update(self) -> None: self.args.extend() self.jh.__init__(**self.args.__dict__) diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 82278bb..ca2d33c 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -169,6 +169,7 @@ def transient_temperature( self, time: floatArray = np.array([]), T0: Optional[float] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -194,49 +195,25 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError("The length of the time array must be at least 2.") - # get initial temperature - if T0 is None: - T0 = ( - self.args.Ta - if isinstance(self.args.Ta, numbers.Number) - else self.args.Ta[0] - ) - - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) - - # Two dicts, one (dc) with static quantities (with all elements of size n), the other (de) - # with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # save args + args = self.args.__dict__.copy() + + # initial conditions T = np.zeros((N, n)) + if T0 is None: + T0 = self.args.Ta T[0, :] = T0 - # main time loop - for i in range(1, len(time)): - for k, v in de.items(): + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): self.args[k] = v[i, :] self.update() T[i, :] = ( @@ -244,15 +221,15 @@ def transient_temperature( ) # save results - dr = dict(time=time, T=T) + dr = {Solver_.Names.time: time, Solver_.Names.temp: T} - # manage return dict 2 : powers + # add power to return dict if needed if return_power: for c in Solver_.Names.powers(): dr[c] = np.zeros_like(T) for i in range(N): - for k in de.keys(): - self.args[k] = de[k][i, :] + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() dr[Solver_.Names.pjle][i, :] = self.jh.value(T[i, :]) dr[Solver_.Names.psol][i, :] = self.sh.value(T[i, :]) @@ -260,13 +237,14 @@ def transient_temperature( dr[Solver_.Names.prad][i, :] = self.rc.value(T[i, :]) dr[Solver_.Names.ppre][i, :] = self.pc.value(T[i, :]) - # squeeze return values if n is 1 + # squeeze values in return dict (if n is 1) if n == 1: keys = list(dr.keys()) keys.remove(Solver_.Names.time) for k in keys: dr[k] = dr[k][:, 0] - return dr - + # restore args + self.args = Args(args) + return dr diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 22fbc28..023a547 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -12,7 +12,7 @@ from thermohl import floatArrayLike, floatArray, strListLike, intArray from thermohl.power import PowerTerm -from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args, _set_dates, reshape +from thermohl.solver.base import Solver as Solver_, _DEFPARAM as DP, Args from thermohl.solver.slv1t import Solver1T from thermohl.utils import quasi_newton_2d @@ -501,6 +501,7 @@ def transient_temperature( time: floatArray = np.array([]), Ts0: Optional[floatArrayLike] = None, Tc0: Optional[floatArrayLike] = None, + dynamic: dict = None, return_power: bool = False, ) -> Dict[str, Any]: """ @@ -531,51 +532,31 @@ def transient_temperature( # get sizes (n for input dict entries, N for time) n = self._min_shape()[0] N = len(time) - if N < 2: - raise ValueError() - - # get initial temperature - Ts0 = Ts0 if Ts0 is not None else self.args.Ta - Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 - # get month, day and hours - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) - # Two dicts, one (dc) with static quantities (with all elements of size - # n), the other (de) with time-changing quantities (with all elements of - # size N*n); uk is a list of keys that are in dc but not in de. - de = dict( - month=month, - day=day, - hour=hour, - I=reshape(self.args.I, N, n), - Ta=reshape(self.args.Ta, N, n), - wa=reshape(self.args.wa, N, n), - ws=reshape(self.args.ws, N, n), - Pa=reshape(self.args.Pa, N, n), - rh=reshape(self.args.rh, N, n), - pr=reshape(self.args.pr, N, n), - ) - del (month, day, hour) + # save args + args = self.args.__dict__.copy() # shortcuts for time-loop c1, c2 = self._morgan_transient() imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) - # main time loop - for i in range(1, len(time)): - for k in de.keys(): - self.args[k] = de[k][i, :] + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): + self.args[k] = v[i, :] self.update() bal = self.balance(ts[i - 1, :], tc[i - 1, :]) ta[i, :] = ta[i - 1, :] + (time[i] - time[i - 1]) * bal * imc @@ -583,4 +564,10 @@ def transient_temperature( tc[i, :] = ta[i, :] + c2 * mrg ts[i, :] = tc[i, :] - mrg - return self._transient_temperature_results(time, ts, ta, tc, return_power, n) + # get results + dr = self._transient_temperature_results(time, ts, ta, tc, return_power, n) + + # restore args + self.args = Args(args) + + return dr From bfcb50f3886d409db805a598f415a2f33482758d Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:30:25 +0200 Subject: [PATCH 42/53] Updated tests. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv3t_legacy.py | 6 +- test/functional_test/functional_test.py | 110 +++++++++++++++++++++++- test/unit/solver/test_slv1t.py | 14 +-- 3 files changed, 118 insertions(+), 12 deletions(-) diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 53597f0..b4171a0 100644 --- a/src/thermohl/solver/slv3t_legacy.py +++ b/src/thermohl/solver/slv3t_legacy.py @@ -178,17 +178,19 @@ def transient_temperature_legacy( # shortcuts for time-loop imc = 1.0 / (self.args.m * self.args.c) - # init + # initial conditions ts = np.zeros((N, n)) ta = np.zeros((N, n)) tc = np.zeros((N, n)) dT = np.zeros((N, n)) - + Ts0 = Ts0 if Ts0 is not None else self.args.Ta + Tc0 = Tc0 if Tc0 is not None else 1.0 + Ts0 ts[0, :] = Ts0 tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) dT[0, :] = tc[0, :] - ts[0, :] + # time loop for i in range(1, len(time)): bal = self.balance(ts[i - 1, :], tc[i - 1, :]) dti = time[i] - time[i - 1] diff --git a/test/functional_test/functional_test.py b/test/functional_test/functional_test.py index 19e927b..b9c164a 100644 --- a/test/functional_test/functional_test.py +++ b/test/functional_test/functional_test.py @@ -10,6 +10,7 @@ import numpy as np import pandas as pd +import thermohl.utils import yaml from typing import List, Dict @@ -67,13 +68,37 @@ def test_steady_temperature(): for d in scenario("temperature", "steady"): for _, e in d.items(): s = cner(scn2dict(e), heateq="3tl") - r = s.steady_temperature() + r = s.steady_temperature(return_power=False) assert np.allclose(r["t_surf"], e["T_surf"], atol=0.05) assert np.allclose(r["t_avg"], e["T_mean"], atol=0.05) assert np.allclose(r["t_core"], e["T_heart"], atol=0.05) +def test_steady_temperature_batch(): + din = [] + dref = [] + for d in scenario("temperature", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("T_surf", "T_mean", "T_heart") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + dref.rename( + columns={"T_surf": "t_surf", "T_mean": "t_avg", "T_heart": "t_core"}, + inplace=True, + ) + + s = cner(din, heateq="3tl") + r = s.steady_temperature(return_power=False) + assert np.allclose(r.values, dref.values, atol=0.05) + + def test_steady_ampacity(): for d in scenario("ampacity", "steady"): for _, e in d.items(): @@ -83,9 +108,29 @@ def test_steady_ampacity(): assert np.allclose(r["I"], e["I_max"], atol=0.05) +def test_steady_ampacity_batch(): + din = [] + dref = [] + for d in scenario("ampacity", "steady"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + dref.append( + pd.DataFrame( + {k: e[k] for k in ("Tmax_cable", "I_max") if k in e}, + index=(0,), + ) + ) + din = pd.concat(din).reset_index(drop=True) + dref = pd.concat(dref).reset_index(drop=True) + + s = cner(din, heateq="3tl") + r = s.steady_intensity(T=dref["Tmax_cable"], return_power=False) + assert np.allclose(r["I"].values, dref["I_max"].values, atol=0.05) + + def test_transient_temperature(): - atol = 0.5 + atol = 0.1 # this is hard-coded, maybe it should be put in the yaml file ... tau = 600.0 @@ -109,7 +154,7 @@ def test_transient_temperature(): rf = s.steady_temperature(Tsg=e["T_mean_final"], Tcg=e["T_mean_final"]) # time - time = np.arange(0.0, 1800.0, dt) + time = np.arange(0.0, 1800.0 + dt, dt) # transient temperature (linearized) rl = s.transient_temperature_legacy( @@ -128,3 +173,62 @@ def test_transient_temperature(): expected_temp = np.array(list(e[k1].values())) estimated_temp = np.interp(expected_time, rl["time"], rl[k2]) assert np.allclose(expected_temp, estimated_temp, atol=atol) + + +def test_transient_temperature_batch(): + + atol = 0.1 + + # this is hard-coded, maybe it should be put in the yaml file ... + tau = 600.0 + dt = 10.0 + minute = 60 + + din = [] + I0 = [] + dref = {"time": [], "t_surf": [], "t_avg": [], "t_core": [], "t_avg_final": []} + for d in scenario("temperature", "transient"): + for _, e in d.items(): + din.append(pd.DataFrame(scn2dict(e), index=(0,))) + I0.append(e["I0_cable"]) + for k1, k2 in zip( + ["T_surf_transient", "T_mean_transient", "T_heart_transient"], + ["t_surf", "t_avg", "t_core"], + ): + dref[k2].append(np.array(list(e[k1].values()))) + dref["t_avg_final"].append(e["T_mean_final"]) + din = pd.concat(din).reset_index(drop=True) + I0 = np.array(I0) + dref["time"] = np.array(list(e["T_surf_transient"].keys())) * minute * 1.0 + for k in ["t_surf", "t_avg", "t_core"]: + dref[k] = np.array(dref[k]).T + dref["t_avg_final"] = np.array(dref["t_avg_final"]) + + # solver + s = cner(din, heateq="3tl") + + # initial steady state + s.args["I"] = I0 + s.update() + ri = s.steady_temperature(return_power=False) + + # final steady state + s.args["I"] = din["I"] + s.update() + rf = s.steady_temperature(Tsg=dref["t_avg_final"], Tcg=dref["t_avg_final"]) + + # check final temp + assert np.allclose(dref["t_avg_final"], rf["t_avg"].values, atol=atol) + + # time + time = np.arange(0.0, 1800.0 + dt, dt) + + # transient temperature (linearized) + rl = s.transient_temperature_legacy( + time=time, Ts0=ri["t_surf"], Tc0=ri["t_core"], tau=tau + ) + + # filter expected times + ix = np.in1d(rl["time"], dref["time"]).nonzero()[0] + for k in ["t_surf", "t_avg", "t_core"]: + assert np.allclose(dref[k], rl[k][ix, :], atol=atol) diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 3cbe419..f9626d9 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -87,9 +87,9 @@ def test_transient_temperature_default(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) def test_transient_temperature_with_initial_temp(solver): @@ -100,10 +100,10 @@ def test_transient_temperature_with_initial_temp(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) - assert result["T"][0] == T0 + assert len(result["t"]) == len(time) + assert result["t"][0] == T0 def test_transient_temperature_with_error(solver): @@ -113,9 +113,9 @@ def test_transient_temperature_with_error(solver): assert isinstance(result, dict) assert "time" in result - assert "T" in result + assert "t" in result assert len(result["time"]) == len(time) - assert len(result["T"]) == len(time) + assert len(result["t"]) == len(time) assert Solver1T.Names.pjle in result assert Solver1T.Names.psol in result assert Solver1T.Names.pcnv in result From b80e1b43c50f2b3dc186a2a5dd1a09dcadf49eed Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:33:28 +0200 Subject: [PATCH 43/53] Updated tests. Signed-off-by: Emmanuel Cieren --- test/test_shape.py | 21 +++++++++++++-------- test/test_solver3t.py | 1 + 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/test/test_shape.py b/test/test_shape.py index 20be2f4..243803a 100644 --- a/test/test_shape.py +++ b/test/test_shape.py @@ -29,14 +29,15 @@ def _ampargs(s: solver.Solver, t: pd.DataFrame): return a -def _traargs(s: solver.Solver, ds: pd.DataFrame, t, I): +def _traargs(s: solver.Solver, ds: pd.DataFrame, t: np.ndarray, dyn: dict): if isinstance(s, solver.Solver1T): - a = dict(time=t, T0=ds[solver.Solver.Names.temp].values) + a = dict(time=t, T0=ds[solver.Solver.Names.temp].values, dynamic=dyn) elif isinstance(s, solver.Solver3T): a = dict( time=t, Ts0=ds[solver.Solver.Names.tsurf].values, Tc0=ds[solver.Solver.Names.tcore].values, + dynamic=dyn, ) else: raise NotImplementedError @@ -104,13 +105,11 @@ def test_steady_1d_mix(): assert len(i) == n -def test_transient_0(): +def test_transient_no_dyn(): for s in _solvers(): t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) - ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, None) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) @@ -131,10 +130,16 @@ def test_transient_1(): s.update() t = np.linspace(0, 3600, 361) - I = 199 * np.ones_like(t) + d = { + "I": 199 * np.ones_like(t), + "Ta": 33.0, + "Pa": np.array(1.013e05), + "ws": 2.0 * np.ones(n), + "wa": np.ones((len(t), n)), + } ds = s.steady_temperature() - a = _traargs(s, ds, t, I) + a = _traargs(s, ds, t, d) r = s.transient_temperature(**a) assert len(r.pop("time")) == len(t) diff --git a/test/test_solver3t.py b/test/test_solver3t.py index 5dd4b73..fe2bc9f 100644 --- a/test/test_solver3t.py +++ b/test/test_solver3t.py @@ -86,6 +86,7 @@ def test_consistency(): target=t, return_err=True, return_power=True, + return_temp=True, tol=1.0e-09, maxiter=64, ) From e3110af81a8366faad0dc53d7cce34910c6f2492 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 08:53:12 +0200 Subject: [PATCH 44/53] Changed max_len/exten/compress methods in Solver class; updated tests Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 35 ++--------------------------------- src/thermohl/solver/slv1t.py | 1 - 2 files changed, 2 insertions(+), 34 deletions(-) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 7eb4d39..3e1d4d6 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -122,7 +122,7 @@ def shape(self) -> Tuple[int, ...]: Compute the maximum effective shape of values in current instance. Members of Args can be float of arrays. If arrays, they must be - one-diemnsional. Float and 1d array can coexist, but all arrays should + one-dimensional. Float and 1d array can coexist, but all arrays should have the same shape/size. This method iterates over all keys in the instance's __dict__ and @@ -260,37 +260,6 @@ def _min_shape(self) -> Tuple[int, ...]: shape = (1,) return shape - def _transient_process_dynamic( - self, time: np.ndarray, n: int, dynamic: dict = None - ): - """Code factorization for transient temperature computations. - - This methods prepare a dict with dynamic values to use in the - compute time loop. - """ - if len(time) < 2: - raise ValueError("The length of the time array must be at least 2.") - - # get month, day and hours in range with time - month, day, hour = _set_dates( - self.args.month, self.args.day, self.args.hour, time, n - ) - - # put dynamic values in a separate dict which will be used - # through the time loop - dynamic_ = { - "month": month, - "day": day, - "hour": hour, - } - - if dynamic is None: - dynamic = {} - for k, v in dynamic.items(): - dynamic_[k] = reshape(v, len(time), n) - - return dynamic_ - def update(self) -> None: self.args.extend() self.jh.__init__(**self.args.__dict__) @@ -327,7 +296,7 @@ def reshape(input_var: numberArrayLike, nb_row: int, nb_columns: int) -> numberA Reshape the input array to the specified dimensions (nr, nc) if possible. Args: - input (numberArrayLike): Variable to be reshaped. + input_var (numberArrayLike): Variable to be reshaped. nb_row (int): Desired number of rows for the reshaped array. nb_columns (int): Desired number of columns for the reshaped array. diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index ca2d33c..594d05d 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -14,7 +14,6 @@ from thermohl import floatArrayLike, floatArray from thermohl.solver.base import Solver as Solver_, Args from thermohl.solver.base import _DEFPARAM as DP -from thermohl.solver.base import _set_dates, reshape from thermohl.utils import bisect_v From 98459deb54cd0af2690d577515a158823dcf4214 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Tue, 1 Apr 2025 09:29:59 +0200 Subject: [PATCH 45/53] Added dynamic parameter to transient temperature methods; factorized some code. Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/base.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 3e1d4d6..2ce8916 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -260,6 +260,37 @@ def _min_shape(self) -> Tuple[int, ...]: shape = (1,) return shape + def _transient_process_dynamic( + self, time: np.ndarray, n: int, dynamic: dict = None + ): + """Code factorization for transient temperature computations. + + This methods prepare a dict with dynamic values to use in the + compute time loop. + """ + if len(time) < 2: + raise ValueError("The length of the time array must be at least 2.") + + # get month, day and hours in range with time + month, day, hour = _set_dates( + self.args.month, self.args.day, self.args.hour, time, n + ) + + # put dynamic values in a separate dict which will be used + # through the time loop + dynamic_ = { + "month": month, + "day": day, + "hour": hour, + } + + if dynamic is None: + dynamic = {} + for k, v in dynamic.items(): + dynamic_[k] = reshape(v, len(time), n) + + return dynamic_ + def update(self) -> None: self.args.extend() self.jh.__init__(**self.args.__dict__) From 0f673b6c23bab81aa90b2d2cb77d5a5490e196a7 Mon Sep 17 00:00:00 2001 From: ecieren Date: Wed, 2 Apr 2025 09:56:55 +0200 Subject: [PATCH 46/53] updated env name Signed-off-by: ecieren --- examples/create_env | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/create_env b/examples/create_env index 68bf9be..f36ab6a 100755 --- a/examples/create_env +++ b/examples/create_env @@ -5,7 +5,7 @@ ENV_DIR=/home/${USER}/ENV mkdir -p $ENV_DIR # -- current env with link -VENV=${ENV_DIR}/thermohl +VENV=${ENV_DIR}/ecieren-thermohl VLNK=.venv ln -sf $VENV $VLNK python3 -m venv $VENV From 1aa55f02a4fe4a639faa18fc23cefb94821bdf7b Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren <129041266+ecieren@users.noreply.github.com> Date: Wed, 2 Apr 2025 10:41:34 +0200 Subject: [PATCH 47/53] Update sonar-project.properties Signed-off-by: Emmanuel Cieren <129041266+ecieren@users.noreply.github.com> --- sonar-project.properties | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sonar-project.properties b/sonar-project.properties index c1ddfd3..eaa7379 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -6,11 +6,11 @@ # SPDX-License-Identifier: MPL-2.0 sonar.host.url=https://sonarcloud.io -sonar.projectKey=ecieren_thermohl -sonar.organization=ecieren +sonar.organization=phlowers +sonar.projectKey=phlowers_thermohl sonar.sources=src/ sonar.tests=test/ sonar.python.version=3.9,3.10,3.11 -sonar.python.coverage.reportPaths=coverage.xml \ No newline at end of file +sonar.python.coverage.reportPaths=coverage.xml From 4360bc9edfe0268ded9ce53f36fd49cc1a753ef9 Mon Sep 17 00:00:00 2001 From: ecieren Date: Mon, 5 May 2025 16:58:26 +0200 Subject: [PATCH 48/53] Formatting (including imports) Signed-off-by: ecieren --- src/thermohl/power/__init__.py | 3 +-- src/thermohl/power/cigre/__init__.py | 2 +- src/thermohl/power/cner/__init__.py | 2 +- src/thermohl/power/ieee/__init__.py | 1 - src/thermohl/power/ieee/convective_cooling.py | 1 - src/thermohl/power/ieee/solar_heating.py | 1 - src/thermohl/power/olla/__init__.py | 2 +- src/thermohl/power/olla/solar_heating.py | 2 +- src/thermohl/solver/__init__.py | 1 - src/thermohl/solver/base.py | 1 - src/thermohl/solver/slv1t.py | 1 - src/thermohl/sun.py | 1 + src/thermohl/thermodynamics.py | 1 + src/thermohl/uncertainties.py | 1 + test/functional_test/functional_test.py | 3 +-- test/test_ieee.py | 2 +- test/unit/power/cigre/test_power_cigre_air.py | 2 +- test/unit/power/cigre/test_power_cigre_convective_cooling.py | 3 +-- test/unit/power/cner/test_power_cner_convective_cooling.py | 3 +-- test/unit/power/cner/test_power_cner_joule_heating.py | 2 +- test/unit/power/ieee/test_power_ieee_joule_heating.py | 2 +- test/unit/power/ieee/test_power_ieee_radiative_cooling.py | 4 ++-- test/unit/power/olla/test_power_olla_air.py | 1 + test/unit/power/test_convective_cooling.py | 2 +- test/unit/power/test_radiative_cooling.py | 2 +- test/unit/solver/test_slv1t.py | 3 ++- 26 files changed, 22 insertions(+), 27 deletions(-) diff --git a/src/thermohl/power/__init__.py b/src/thermohl/power/__init__.py index 6d4c60b..0b8cf04 100644 --- a/src/thermohl/power/__init__.py +++ b/src/thermohl/power/__init__.py @@ -5,7 +5,6 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -from .radiative_cooling import RadiativeCoolingBase - from .power_term import PowerTerm +from .radiative_cooling import RadiativeCoolingBase from .solar_heating import _SRad, SolarHeatingBase diff --git a/src/thermohl/power/cigre/__init__.py b/src/thermohl/power/cigre/__init__.py index b957cb6..2c4d597 100644 --- a/src/thermohl/power/cigre/__init__.py +++ b/src/thermohl/power/cigre/__init__.py @@ -12,7 +12,7 @@ """ from .air import Air -from .solar_heating import SolarHeating from .convective_cooling import ConvectiveCooling from .joule_heating import JouleHeating from .radiative_cooling import RadiativeCooling +from .solar_heating import SolarHeating diff --git a/src/thermohl/power/cner/__init__.py b/src/thermohl/power/cner/__init__.py index 10ef691..2217e32 100644 --- a/src/thermohl/power/cner/__init__.py +++ b/src/thermohl/power/cner/__init__.py @@ -11,7 +11,7 @@ """ from .air import Air -from .solar_heating import SolarHeating from .convective_cooling import ConvectiveCooling from .joule_heating import JouleHeating from .radiative_cooling import RadiativeCooling +from .solar_heating import SolarHeating diff --git a/src/thermohl/power/ieee/__init__.py b/src/thermohl/power/ieee/__init__.py index ee68dbd..9feef25 100644 --- a/src/thermohl/power/ieee/__init__.py +++ b/src/thermohl/power/ieee/__init__.py @@ -13,7 +13,6 @@ from .air import Air from .convective_cooling import ConvectiveCooling -from ..convective_cooling import ConvectiveCoolingBase from .joule_heating import JouleHeating from .radiative_cooling import RadiativeCooling from .solar_heating import SolarHeating diff --git a/src/thermohl/power/ieee/convective_cooling.py b/src/thermohl/power/ieee/convective_cooling.py index 13a326b..e15a49a 100644 --- a/src/thermohl/power/ieee/convective_cooling.py +++ b/src/thermohl/power/ieee/convective_cooling.py @@ -7,7 +7,6 @@ from typing import Any - from thermohl import floatArrayLike from thermohl.power.convective_cooling import ConvectiveCoolingBase from thermohl.power.ieee import Air diff --git a/src/thermohl/power/ieee/solar_heating.py b/src/thermohl/power/ieee/solar_heating.py index 0145670..c41ff7b 100644 --- a/src/thermohl/power/ieee/solar_heating.py +++ b/src/thermohl/power/ieee/solar_heating.py @@ -7,7 +7,6 @@ from typing import Optional, Any - from thermohl import floatArrayLike, intArrayLike from thermohl.power import _SRad, SolarHeatingBase diff --git a/src/thermohl/power/olla/__init__.py b/src/thermohl/power/olla/__init__.py index 64431db..f5d1a28 100644 --- a/src/thermohl/power/olla/__init__.py +++ b/src/thermohl/power/olla/__init__.py @@ -11,7 +11,7 @@ """ from .air import Air -from .solar_heating import SolarHeating from .convective_cooling import ConvectiveCooling from .joule_heating import JouleHeating from .radiative_cooling import RadiativeCooling +from .solar_heating import SolarHeating diff --git a/src/thermohl/power/olla/solar_heating.py b/src/thermohl/power/olla/solar_heating.py index 99bc5a7..e03c444 100644 --- a/src/thermohl/power/olla/solar_heating.py +++ b/src/thermohl/power/olla/solar_heating.py @@ -7,8 +7,8 @@ from typing import Optional, Any -from thermohl.power import ieee from thermohl import floatArrayLike, intArrayLike +from thermohl.power import ieee class SolarHeating(ieee.SolarHeating): diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index 36d2f30..73ec51a 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -13,7 +13,6 @@ from thermohl.power import cner as _cner from thermohl.power import ieee as _ieee from thermohl.power import olla as _olla - from thermohl.solver.base import Args, Solver from thermohl.solver.slv1d import Solver1D from thermohl.solver.slv1t import Solver1T diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index 2ce8916..964cb85 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -13,7 +13,6 @@ import numpy as np import pandas as pd -from numpy import ndarray from thermohl import ( floatArrayLike, diff --git a/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 594d05d..7777ef2 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -5,7 +5,6 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import numbers from typing import Dict, Any, Optional import numpy as np diff --git a/src/thermohl/sun.py b/src/thermohl/sun.py index 3b99fde..4feaec9 100644 --- a/src/thermohl/sun.py +++ b/src/thermohl/sun.py @@ -13,6 +13,7 @@ """ import numpy as np + from thermohl import floatArrayLike, intArrayLike diff --git a/src/thermohl/thermodynamics.py b/src/thermohl/thermodynamics.py index 201d80d..c191bd9 100644 --- a/src/thermohl/thermodynamics.py +++ b/src/thermohl/thermodynamics.py @@ -8,6 +8,7 @@ """Thermodynamics quantities. Values from Wikipedia unless specified.""" import numpy as np + from thermohl import floatArrayLike # Standard temperature and pressure from EPA and NIST; _T0 in K, _p0 in Pa diff --git a/src/thermohl/uncertainties.py b/src/thermohl/uncertainties.py index eb3db4d..e9ae113 100644 --- a/src/thermohl/uncertainties.py +++ b/src/thermohl/uncertainties.py @@ -15,6 +15,7 @@ import scipy.stats from scipy.stats import circmean from scipy.stats import circstd + from thermohl import distributions, frozen_dist from thermohl import solver from thermohl import utils diff --git a/test/functional_test/functional_test.py b/test/functional_test/functional_test.py index b9c164a..da0aeaa 100644 --- a/test/functional_test/functional_test.py +++ b/test/functional_test/functional_test.py @@ -7,12 +7,11 @@ import datetime import os.path +from typing import List, Dict import numpy as np import pandas as pd -import thermohl.utils import yaml -from typing import List, Dict from thermohl.solver import cner diff --git a/test/test_ieee.py b/test/test_ieee.py index 5637f31..c6d7e13 100644 --- a/test/test_ieee.py +++ b/test/test_ieee.py @@ -7,8 +7,8 @@ import numpy as np -from thermohl.power import ieee from thermohl import solver +from thermohl.power import ieee def test_compare_powers(): diff --git a/test/unit/power/cigre/test_power_cigre_air.py b/test/unit/power/cigre/test_power_cigre_air.py index 8543121..8f5d890 100644 --- a/test/unit/power/cigre/test_power_cigre_air.py +++ b/test/unit/power/cigre/test_power_cigre_air.py @@ -5,8 +5,8 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl.power.cigre import Air 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 bc51875..9ece8d1 100644 --- a/test/unit/power/cigre/test_power_cigre_convective_cooling.py +++ b/test/unit/power/cigre/test_power_cigre_convective_cooling.py @@ -8,9 +8,8 @@ import numpy as np import pytest -from thermohl.power.cigre import ConvectiveCooling from thermohl.power.cigre import Air - +from thermohl.power.cigre import ConvectiveCooling conv_cool_instances = [ ( diff --git a/test/unit/power/cner/test_power_cner_convective_cooling.py b/test/unit/power/cner/test_power_cner_convective_cooling.py index e1fecd6..0eb7033 100644 --- a/test/unit/power/cner/test_power_cner_convective_cooling.py +++ b/test/unit/power/cner/test_power_cner_convective_cooling.py @@ -5,12 +5,11 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl.power.cner import ConvectiveCooling - conv_cool_instances = [ ConvectiveCooling( alt=np.array([100.0]), diff --git a/test/unit/power/cner/test_power_cner_joule_heating.py b/test/unit/power/cner/test_power_cner_joule_heating.py index cd702d7..1985049 100644 --- a/test/unit/power/cner/test_power_cner_joule_heating.py +++ b/test/unit/power/cner/test_power_cner_joule_heating.py @@ -5,8 +5,8 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl.power.cner import JouleHeating diff --git a/test/unit/power/ieee/test_power_ieee_joule_heating.py b/test/unit/power/ieee/test_power_ieee_joule_heating.py index 069e53b..fa85309 100644 --- a/test/unit/power/ieee/test_power_ieee_joule_heating.py +++ b/test/unit/power/ieee/test_power_ieee_joule_heating.py @@ -5,8 +5,8 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl.power.ieee import JouleHeating diff --git a/test/unit/power/ieee/test_power_ieee_radiative_cooling.py b/test/unit/power/ieee/test_power_ieee_radiative_cooling.py index 840f39a..b2718ca 100644 --- a/test/unit/power/ieee/test_power_ieee_radiative_cooling.py +++ b/test/unit/power/ieee/test_power_ieee_radiative_cooling.py @@ -5,11 +5,11 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest -from thermohl.power.ieee import RadiativeCooling from thermohl import solver +from thermohl.power.ieee import RadiativeCooling def set_default_values_scalar(): diff --git a/test/unit/power/olla/test_power_olla_air.py b/test/unit/power/olla/test_power_olla_air.py index 95a10a9..381cbe3 100644 --- a/test/unit/power/olla/test_power_olla_air.py +++ b/test/unit/power/olla/test_power_olla_air.py @@ -6,6 +6,7 @@ # SPDX-License-Identifier: MPL-2.0 import numpy as np + from thermohl.power.olla import Air diff --git a/test/unit/power/test_convective_cooling.py b/test/unit/power/test_convective_cooling.py index c104559..864205f 100644 --- a/test/unit/power/test_convective_cooling.py +++ b/test/unit/power/test_convective_cooling.py @@ -5,8 +5,8 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl import solver from thermohl.power.ieee import ConvectiveCooling diff --git a/test/unit/power/test_radiative_cooling.py b/test/unit/power/test_radiative_cooling.py index baf9ca9..1cfb450 100644 --- a/test/unit/power/test_radiative_cooling.py +++ b/test/unit/power/test_radiative_cooling.py @@ -5,8 +5,8 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np +import pytest from thermohl.power.radiative_cooling import RadiativeCoolingBase diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index f9626d9..b49ea9c 100644 --- a/test/unit/solver/test_slv1t.py +++ b/test/unit/solver/test_slv1t.py @@ -5,9 +5,10 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -import pytest import numpy as np import pandas as pd +import pytest + from thermohl.solver.slv1t import Solver1T From 285b9d7da9f66a961766b3328749a90e59733ad8 Mon Sep 17 00:00:00 2001 From: ecieren Date: Mon, 19 May 2025 11:14:54 +0200 Subject: [PATCH 49/53] fix test Signed-off-by: ecieren --- test/test_power_rte_integration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_power_rte_integration.py b/test/test_power_rte_integration.py index a9cbcff..6242e05 100644 --- a/test/test_power_rte_integration.py +++ b/test/test_power_rte_integration.py @@ -220,6 +220,6 @@ def test_solar_heating(): alpha = 0.9 * ones p = np.array([34.9, 21.9357, 13.95, 21.9357, 21.9357]) - s = rte.SolarHeating(lat, azm, month, day, hour, D, alpha, srad=None) + s = rte.SolarHeating(lat, azm, month, day, hour, D, alpha) assert np.allclose(p, s.value(ones), 0.1) From b4501b35a5a51ba31b90aba290cf3f435d8aa6df Mon Sep 17 00:00:00 2001 From: ecieren Date: Mon, 19 May 2025 11:34:23 +0200 Subject: [PATCH 50/53] Formatting (including imports) Signed-off-by: ecieren --- src/thermohl/power/olla/__init__.py | 3 +-- src/thermohl/power/rte/__init__.py | 3 +-- src/thermohl/power/rte/convective_cooling.py | 2 +- src/thermohl/solver/__init__.py | 6 ++++-- test/functional_test/functional_test.py | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/thermohl/power/olla/__init__.py b/src/thermohl/power/olla/__init__.py index 86684e4..6b94a0c 100644 --- a/src/thermohl/power/olla/__init__.py +++ b/src/thermohl/power/olla/__init__.py @@ -5,8 +5,7 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -"""Power terms implementation using RTE's olla project choices. -""" +"""Power terms implementation using RTE's olla project choices.""" from .air import Air from .convective_cooling import ConvectiveCooling diff --git a/src/thermohl/power/rte/__init__.py b/src/thermohl/power/rte/__init__.py index 371980a..90da7b7 100644 --- a/src/thermohl/power/rte/__init__.py +++ b/src/thermohl/power/rte/__init__.py @@ -5,8 +5,7 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. # SPDX-License-Identifier: MPL-2.0 -"""Power terms implementation matching rte's Excel sheet. -""" +"""Power terms implementation matching rte's Excel sheet.""" from .air import Air from .convective_cooling import ConvectiveCooling diff --git a/src/thermohl/power/rte/convective_cooling.py b/src/thermohl/power/rte/convective_cooling.py index 5c366e2..b5103a1 100644 --- a/src/thermohl/power/rte/convective_cooling.py +++ b/src/thermohl/power/rte/convective_cooling.py @@ -10,8 +10,8 @@ import numpy as np from thermohl import floatArrayLike -from thermohl.power.rte import Air from thermohl.power.convective_cooling import ConvectiveCoolingBase +from thermohl.power.rte import Air class ConvectiveCooling(ConvectiveCoolingBase): diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index 484009a..2a6d782 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -10,16 +10,18 @@ from typing import Dict, Any, Optional, Union, Type from thermohl.power import cigre as _cigre -from thermohl.power import rte as _rte from thermohl.power import ieee as _ieee from thermohl.power import olla as _olla +from thermohl.power import rte as _rte from thermohl.solver.base import Args, Solver from thermohl.solver.slv1d import Solver1D from thermohl.solver.slv1t import Solver1T from thermohl.solver.slv3t import Solver3T from thermohl.solver.slv3t_legacy import Solver3TL -concreteSolverType = Union[Type[Solver1T], Type[Solver3T], Type[Solver3TL], Type[Solver1D]] +concreteSolverType = Union[ + Type[Solver1T], Type[Solver3T], Type[Solver3TL], Type[Solver1D] +] def default_values() -> Dict[str, Any]: diff --git a/test/functional_test/functional_test.py b/test/functional_test/functional_test.py index 5f8a33f..ec7b7b2 100644 --- a/test/functional_test/functional_test.py +++ b/test/functional_test/functional_test.py @@ -7,11 +7,11 @@ import datetime import os.path +from typing import List, Dict import numpy as np import pandas as pd import yaml -from typing import List, Dict from thermohl.solver import rte From 8136a035898ab3cbbac6261f74e5351a856dd306 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Mon, 23 Jun 2025 22:02:42 +0200 Subject: [PATCH 51/53] corrections for merge Signed-off-by: Emmanuel Cieren --- src/thermohl/sun.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/thermohl/sun.py b/src/thermohl/sun.py index 2834870..1904ac3 100644 --- a/src/thermohl/sun.py +++ b/src/thermohl/sun.py @@ -17,8 +17,8 @@ from thermohl import floatArrayLike, intArrayLike - -def utc2solar_hour(hour, minute=0., second=0., lon=0.): +def utc2solar_hour(hour: floatArrayLike, minute: floatArrayLike = 0., second: floatArrayLike = 0., + lon: floatArrayLike = 0.): """convert utc hour to solar hour adding the longitude contribution If more than one input are numpy arrays, they should have the same size. @@ -47,7 +47,7 @@ def utc2solar_hour(hour, minute=0., second=0., lon=0.): def hour_angle( - hour: floatArrayLike, minute: floatArrayLike = 0.0, second: floatArrayLike = 0.0 + hour: floatArrayLike, minute: floatArrayLike = 0.0, second: floatArrayLike = 0.0 ) -> floatArrayLike: """Compute hour angle. @@ -72,7 +72,8 @@ def hour_angle( return np.radians(15.0 * (solar_hour - 12.0)) -_csm = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) +# Cumulative days of the year at start of each month +_CSM = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) def solar_declination(month: intArrayLike, day: intArrayLike) -> floatArrayLike: @@ -93,7 +94,7 @@ def solar_declination(month: intArrayLike, day: intArrayLike) -> floatArrayLike: Solar declination in radians. """ - doy = _csm[month - 1] + day + doy = _CSM[month - 1] + day return np.deg2rad(23.46) * np.sin(2.0 * np.pi * (doy + 284) / 365.0) From 5168228b32d409989e193c376f08b2035c087a40 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Mon, 23 Jun 2025 22:03:04 +0200 Subject: [PATCH 52/53] updated 3t model for transient temperature Signed-off-by: Emmanuel Cieren --- src/thermohl/solver/slv3t.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index 023a547..84fdb05 100644 --- a/src/thermohl/solver/slv3t.py +++ b/src/thermohl/solver/slv3t.py @@ -238,6 +238,16 @@ def morgan(self, ts: floatArray, tc: floatArray) -> floatArray: c, _, _, _ = self.mgc return (tc - ts) - c * self.joule(ts, tc) / (2.0 * np.pi * self.args.l) + def tau(self, ts: floatArray, tc: floatArray, dt=1.0e-05) -> floatArrayLike: + """Estimation of a time-constant by linearization of the EDO.""" + db = ( + self.balance(ts + dt, tc) + - self.balance(ts - dt, tc) + + self.balance(ts, tc + dt) + - self.balance(ts, tc - dt) + ) / (2 * dt) + return -(self.args.m * self.args.c) / db + # ========================================================================== def _steady_return_opt( @@ -553,16 +563,23 @@ def transient_temperature( tc[0, :] = Tc0 ta[0, :] = self.average(ts[0, :], tc[0, :]) + tx = tc[0, :] - ts[0, :] + ty = c2 * ts[0, :] + (1 - c2) * tc[0, :] + # time loop for i in range(1, N): for k, v in dynamic_.items(): - self.args[k] = v[i, :] + self.args[k] = v[i - 1, :] self.update() + dt = time[i] - time[i - 1] bal = self.balance(ts[i - 1, :], tc[i - 1, :]) - ta[i, :] = ta[i - 1, :] + (time[i] - time[i - 1]) * bal * imc - mrg = c1 * (self.jh.value(ta[i, :]) - bal) - tc[i, :] = ta[i, :] + c2 * mrg - ts[i, :] = tc[i, :] - mrg + tau = self.tau(ts[i - 1, :], tc[i - 1, :]) + ta[i, :] = ta[i - 1, :] + dt * bal * imc + morgan = c1 * (self.jh.value(ta[i, :]) - bal) + tx = tx + dt * (-tx + morgan) / (tau * 0.3) + ty = ty + dt * (-ty + ta[i - 1, :]) / (tau * 0.02) + tc[i, :] = c2 * tx + ty + ts[i, :] = ty - (1 - c2) * tx # get results dr = self._transient_temperature_results(time, ts, ta, tc, return_power, n) From 221c83579b4d34151a0c5f58f27676434d66cc62 Mon Sep 17 00:00:00 2001 From: Emmanuel Cieren Date: Mon, 23 Jun 2025 22:09:58 +0200 Subject: [PATCH 53/53] black formater Signed-off-by: Emmanuel Cieren --- src/thermohl/sun.py | 13 ++++++++----- test/unit/test_sun.py | 27 ++++++++++++++++----------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/thermohl/sun.py b/src/thermohl/sun.py index 1904ac3..c54c29e 100644 --- a/src/thermohl/sun.py +++ b/src/thermohl/sun.py @@ -17,8 +17,12 @@ from thermohl import floatArrayLike, intArrayLike -def utc2solar_hour(hour: floatArrayLike, minute: floatArrayLike = 0., second: floatArrayLike = 0., - lon: floatArrayLike = 0.): +def utc2solar_hour( + hour: floatArrayLike, + minute: floatArrayLike = 0.0, + second: floatArrayLike = 0.0, + lon: floatArrayLike = 0.0, +): """convert utc hour to solar hour adding the longitude contribution If more than one input are numpy arrays, they should have the same size. @@ -41,13 +45,12 @@ def utc2solar_hour(hour: floatArrayLike, minute: floatArrayLike = 0., second: fl """ # add 4 min (1/15 of an hour) for every degree of east longitude - solar_hour = hour % 24 + minute / 60. + second / 3600. + np.rad2deg(lon) / 15. + solar_hour = hour % 24 + minute / 60.0 + second / 3600.0 + np.rad2deg(lon) / 15.0 return solar_hour - def hour_angle( - hour: floatArrayLike, minute: floatArrayLike = 0.0, second: floatArrayLike = 0.0 + hour: floatArrayLike, minute: floatArrayLike = 0.0, second: floatArrayLike = 0.0 ) -> floatArrayLike: """Compute hour angle. diff --git a/test/unit/test_sun.py b/test/unit/test_sun.py index e8d9f17..f5901e9 100644 --- a/test/unit/test_sun.py +++ b/test/unit/test_sun.py @@ -17,9 +17,10 @@ def test_scalar_input(): lon = np.deg2rad(45) # 45 degrees east result = utc2solar_hour(hour, minute, second, lon) - expected_result = 12 + 30 / 60. + 45 / 3600. + 45 / 15. + expected_result = 12 + 30 / 60.0 + 45 / 3600.0 + 45 / 15.0 assert np.isclose(result, expected_result) + def test_array_input(): # Test with numpy array inputs hours = np.array([12, 15, 18]) @@ -28,13 +29,16 @@ def test_array_input(): lons = np.deg2rad(np.array([45, 90, 135])) # 45, 90, and 135 degrees east result = utc2solar_hour(hours, minutes, seconds, lons) - expected_result = np.array([ - 12 + 30 / 60. + 45 / 3600. + 45 / 15., - 15 + 45 / 60. + 30 / 3600. + 90 / 15., - 18 + 0 / 60. + 15 / 3600. + 135 / 15. - ]) + expected_result = np.array( + [ + 12 + 30 / 60.0 + 45 / 3600.0 + 45 / 15.0, + 15 + 45 / 60.0 + 30 / 3600.0 + 90 / 15.0, + 18 + 0 / 60.0 + 15 / 3600.0 + 135 / 15.0, + ] + ) np.testing.assert_array_almost_equal(result, expected_result) + def test_edge_cases(): # Test edge cases hour = 23 @@ -43,7 +47,7 @@ def test_edge_cases(): lon = np.deg2rad(180) # 180 degrees east result = utc2solar_hour(hour, minute, second, lon) - expected_result = 23 + 59 / 60. + 59 / 3600. + 180 / 15. + expected_result = 23 + 59 / 60.0 + 59 / 3600.0 + 180 / 15.0 assert np.isclose(result, expected_result) hour = 0 @@ -52,9 +56,10 @@ def test_edge_cases(): lon = np.deg2rad(-180) # 180 degrees west result = utc2solar_hour(hour, minute, second, lon) - expected_result = 0 + 0 / 60. + 0 / 3600. - 180 / 15. + expected_result = 0 + 0 / 60.0 + 0 / 3600.0 - 180 / 15.0 assert np.isclose(result, expected_result) + def test_realistic_cases(): # Test edge cases hour = 23 * np.ones(3) @@ -63,8 +68,8 @@ def test_realistic_cases(): lon = np.deg2rad(np.array([2.33472, 7.75, -4.48])) # Paris, Strasbourg, Brest result = utc2solar_hour(hour, minute, second, lon) - expected_result = 23 + 59 / 60. + 59 / 3600. - + expected_result = 23 + 59 / 60.0 + 59 / 3600.0 + assert result[0] > expected_result assert result[1] > result[0] - assert result[2] < expected_result \ No newline at end of file + assert result[2] < expected_result