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 diff --git a/sonar-project.properties b/sonar-project.properties index 43ce9a2..eaa7379 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -13,4 +13,4 @@ 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 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/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/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 24f00f9..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 @@ -24,7 +23,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 +51,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/__init__.py b/src/thermohl/power/olla/__init__.py index d893624..6b94a0c 100644 --- a/src/thermohl/power/olla/__init__.py +++ b/src/thermohl/power/olla/__init__.py @@ -5,11 +5,10 @@ # 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 .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 1df8d69..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): @@ -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/src/thermohl/power/rte/__init__.py b/src/thermohl/power/rte/__init__.py index f3a231c..d378331 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 .solar_heating import SolarHeating, solar_irradiance 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/power/rte/solar_heating.py b/src/thermohl/power/rte/solar_heating.py index c9024b9..685c1d0 100644 --- a/src/thermohl/power/rte/solar_heating.py +++ b/src/thermohl/power/rte/solar_heating.py @@ -66,7 +66,7 @@ def __init__( hour: floatArrayLike, D: floatArrayLike, alpha: floatArrayLike, - srad: Optional[floatArrayLike] = None, + srad: Optional[floatArrayLike] = float("nan"), **kwargs: Any, ): r"""Build with args. @@ -90,8 +90,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/solar_heating.py b/src/thermohl/power/solar_heating.py index 134b9bb..25a3081 100644 --- a/src/thermohl/power/solar_heating.py +++ b/src/thermohl/power/solar_heating.py @@ -92,11 +92,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) diff --git a/src/thermohl/solver/__init__.py b/src/thermohl/solver/__init__.py index 1c6091d..2a6d782 100644 --- a/src/thermohl/solver/__init__.py +++ b/src/thermohl/solver/__init__.py @@ -9,18 +9,19 @@ from typing import Dict, Any, Optional, Union, Type -from thermohl.power import cigre as cigrep -from thermohl.power import rte as rtep -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 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]: @@ -45,34 +46,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 == "rte": return solver( dic, - rtep.JouleHeating, - rtep.SolarHeating, - rtep.ConvectiveCooling, - rtep.RadiativeCooling, + _rte.JouleHeating, + _rte.SolarHeating, + _rte.ConvectiveCooling, + _rte.RadiativeCooling, ) else: raise ValueError() diff --git a/src/thermohl/solver/base.py b/src/thermohl/solver/base.py index bc13416..964cb85 100644 --- a/src/thermohl/solver/base.py +++ b/src/thermohl/solver/base.py @@ -9,11 +9,10 @@ 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 -from numpy import ndarray from thermohl import ( floatArrayLike, @@ -53,6 +52,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.""" @@ -62,9 +63,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 +74,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) @@ -115,53 +116,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-dimensional. 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 + This method iterates over all keys in the instance's __dict__ and + computes the maximum length of the values associated with those keys. - 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 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 +172,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 +245,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 +253,45 @@ 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 _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_to_max_len() + self.args.extend() self.jh.__init__(**self.args.__dict__) self.sh.__init__(**self.args.__dict__) self.cc.__init__(**self.args.__dict__) @@ -285,12 +321,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_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. @@ -299,24 +335,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 @@ -362,10 +403,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/src/thermohl/solver/slv1t.py b/src/thermohl/solver/slv1t.py index 03f9cfc..7777ef2 100644 --- a/src/thermohl/solver/slv1t.py +++ b/src/thermohl/solver/slv1t.py @@ -5,20 +5,44 @@ # 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 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 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, @@ -56,128 +80,20 @@ 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 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 - def transient_temperature( - self, - time: floatArray = np.array([]), - T0: Optional[float] = None, - return_power: bool = False, - ) -> Dict[str, Any]: - """ - Compute transient-state temperature. - - Parameters - ---------- - time : numpy.ndarray - A 1D array with times (in seconds) when the temperature needs to be - computed. The array must contain increasing values (undefined - behaviour otherwise). - T0 : float - Initial temperature. If set to None, the ambient temperature from - internal dict will be used. The default is None. - return_power : bool, optional - Return power term values. The default is False. - - Returns : Dict[str, Any] - A dictionary with temperature and other results (depending on inputs) - in the keys. - """ - - # get sizes (n for input dict entries, N for time) - n = self.args.max_len() - 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) - - # shortcuts for time-loop - imc = 1.0 / (self.args.m * self.args.c) - - # init - T = np.zeros((N, n)) - T[0, :] = T0 - - # main time loop - for i in range(1, len(time)): - for k, v in de.items(): - self.args[k] = v[i, :] - self.update() - T[i, :] = ( - T[i - 1, :] + (time[i] - time[i - 1]) * self.balance(T[i - 1, :]) * imc - ) - - # save results - dr = dict(time=time, T=T) - - # manage return dict 2 : powers - 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, :] - self.update() - dr[Solver_.Names.pjle][i, :] = self.jh.value(T[i, :]) - dr[Solver_.Names.psol][i, :] = self.sh.value(T[i, :]) - dr[Solver_.Names.pcnv][i, :] = self.cc.value(T[i, :]) - 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 - 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_intensity( self, T: floatArrayLike = np.array([]), @@ -222,7 +138,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_) @@ -243,15 +159,90 @@ def fun(i: floatArray) -> floatArrayLike: # format output df = pd.DataFrame(data=A, columns=[Solver_.Names.transit]) + df = self._steady_return_opt(return_err, return_power, T_, err, df) - if return_err: - df[Solver_.Names.err] = err + return df + + def transient_temperature( + self, + time: floatArray = np.array([]), + T0: Optional[float] = None, + dynamic: dict = None, + return_power: bool = False, + ) -> Dict[str, Any]: + """ + Compute transient-state temperature. + + Parameters + ---------- + time : numpy.ndarray + A 1D array with times (in seconds) when the temperature needs to be + computed. The array must contain increasing values (undefined + behaviour otherwise). + T0 : float + Initial temperature. If set to None, the ambient temperature from + internal dict will be used. The default is None. + return_power : bool, optional + Return power term values. The default is False. + + Returns : Dict[str, Any] + A dictionary with temperature and other results (depending on inputs) + in the keys. + """ + + # get sizes (n for input dict entries, N for time) + n = self._min_shape()[0] + N = len(time) + + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) + + # shortcuts for time-loop + imc = 1.0 / (self.args.m * self.args.c) + + # 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 + + # time loop + for i in range(1, N): + for k, v in dynamic_.items(): + self.args[k] = v[i, :] + self.update() + T[i, :] = ( + T[i - 1, :] + (time[i] - time[i - 1]) * self.balance(T[i - 1, :]) * imc + ) + + # save results + dr = {Solver_.Names.time: time, Solver_.Names.temp: T} + + # add power to return dict if needed 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) + for c in Solver_.Names.powers(): + dr[c] = np.zeros_like(T) + for i in range(N): + 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, :]) + dr[Solver_.Names.pcnv][i, :] = self.cc.value(T[i, :]) + dr[Solver_.Names.prad][i, :] = self.rc.value(T[i, :]) + dr[Solver_.Names.ppre][i, :] = self.pc.value(T[i, :]) - return df + # 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] + + # restore args + self.args = Args(args) + + return dr diff --git a/src/thermohl/solver/slv3t.py b/src/thermohl/solver/slv3t.py index f038a0b..84fdb05 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 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__( @@ -78,7 +133,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 +152,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__) @@ -183,6 +238,130 @@ 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( + 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_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 + 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, @@ -215,7 +394,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) @@ -240,236 +419,10 @@ 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 - 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, ts, ta, tc, return_power, n): - 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 transient_temperature( - self, - time: floatArray = np.array([]), - Ts0: Optional[floatArrayLike] = None, - Tc0: Optional[floatArrayLike] = None, - return_power: bool = False, - ) -> Dict[str, Any]: - """ - Compute transient-state temperature. - - Parameters - ---------- - time : numpy.ndarray - A 1D array with times (in seconds) when the temperature needs to be - computed. The array must contain increasing values (undefined - behaviour otherwise). - Ts0 : float - Initial surface temperature. If set to None, the ambient temperature from - internal dict will be used. The default is None. - Tc0 : float - Initial core temperature. If set to None, the ambient temperature from - internal dict will be used. The default is None. - return_power : bool, optional - Return power term values. The default is False. - - Returns - ------- - Dict[str, Any] - A dictionary with temperature and other results (depending on inputs) - in the keys. - - """ - # get sizes (n for input dict entries, N for time) - n = self.args.max_len() - 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 - ) - - # 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) - - # shortcuts for time-loop - c1, c2 = self._morgan_transient() - imc = 1.0 / (self.args.m * self.args.c) - - # init - ts = np.zeros((N, n)) - ta = np.zeros((N, n)) - tc = np.zeros((N, n)) - 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, :] - self.update() - 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 - - 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.""" - - max_len = self.args.max_len() - Tmax = T * np.ones(max_len) - target_ = self._check_target(target, self.args.d, max_len) - - # 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 steady_intensity( self, T: floatArrayLike = np.array([]), @@ -539,24 +492,99 @@ 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 + + def transient_temperature( + self, + time: floatArray = np.array([]), + Ts0: Optional[floatArrayLike] = None, + Tc0: Optional[floatArrayLike] = None, + dynamic: dict = None, + return_power: bool = False, + ) -> Dict[str, Any]: + """ + Compute transient-state temperature. + + Parameters + ---------- + time : numpy.ndarray + A 1D array with times (in seconds) when the temperature needs to be + computed. The array must contain increasing values (undefined + behaviour otherwise). + Ts0 : float + Initial surface temperature. If set to None, the ambient temperature from + internal dict will be used. The default is None. + Tc0 : float + Initial core temperature. If set to None, the ambient temperature from + internal dict will be used. The default is None. + return_power : bool, optional + Return power term values. The default is False. + + Returns + ------- + Dict[str, Any] + A dictionary with temperature and other results (depending on inputs) + in the keys. + + """ + # get sizes (n for input dict entries, N for time) + n = self._min_shape()[0] + N = len(time) + + # process dynamic values + dynamic_ = self._transient_process_dynamic(time, n, dynamic) + + # 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) + + # 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, :]) + + 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 - 1, :] + self.update() + dt = time[i] - time[i - 1] + bal = self.balance(ts[i - 1, :], tc[i - 1, :]) + 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) + + # restore args + self.args = Args(args) + + return dr diff --git a/src/thermohl/solver/slv3t_legacy.py b/src/thermohl/solver/slv3t_legacy.py index 2da3a2f..b4171a0 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): @@ -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 @@ -86,14 +86,16 @@ 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]: """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_ = _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([]), @@ -162,7 +166,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() @@ -174,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/src/thermohl/sun.py b/src/thermohl/sun.py index 7f09f9a..c54c29e 100644 --- a/src/thermohl/sun.py +++ b/src/thermohl/sun.py @@ -13,11 +13,16 @@ """ import numpy as np -from thermohl import floatArrayLike, intArrayLike +from thermohl import floatArrayLike, intArrayLike -def utc2solar_hour(hour, minute=0., second=0., lon=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. @@ -40,11 +45,10 @@ def utc2solar_hour(hour, minute=0., second=0., lon=0.): """ # 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 ) -> floatArrayLike: @@ -71,7 +75,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: @@ -92,7 +97,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) 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/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/functional_test/functional_test.py b/test/functional_test/functional_test.py index c899145..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 @@ -67,13 +67,37 @@ def test_steady_temperature(): for d in scenario("temperature", "steady"): for _, e in d.items(): s = rte(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 = rte(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 +107,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 = rte(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 +153,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 +172,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 = rte(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/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/test_power_rte_integration.py b/test/test_power_rte_integration.py index 9e0becb..6242e05 100644 --- a/test/test_power_rte_integration.py +++ b/test/test_power_rte_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 = rte.JouleHeating(**d1) @@ -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) diff --git a/test/test_shape.py b/test/test_shape.py index ff53284..04ae2ca 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 @@ -95,7 +96,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) @@ -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 f4f7848..e47fb47 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, ) 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/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/rte/test_power_rte_convective_cooling.py b/test/unit/power/rte/test_power_rte_convective_cooling.py index f94eaef..bef8657 100644 --- a/test/unit/power/rte/test_power_rte_convective_cooling.py +++ b/test/unit/power/rte/test_power_rte_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.rte import ConvectiveCooling - conv_cool_instances = [ ConvectiveCooling( alt=np.array([100.0]), diff --git a/test/unit/power/rte/test_power_rte_joule_heating.py b/test/unit/power/rte/test_power_rte_joule_heating.py index bb755e1..5750a5b 100644 --- a/test/unit/power/rte/test_power_rte_joule_heating.py +++ b/test/unit/power/rte/test_power_rte_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.rte import JouleHeating 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_base.py b/test/unit/solver/test_base.py index 63eb841..f3822aa 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() - - assert result == 3 + try: + args = Args(dic) + assert False + except ValueError: + pass -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 @@ -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]) @@ -256,7 +267,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 +290,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 +313,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 +339,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 +351,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 diff --git a/test/unit/solver/test_slv1t.py b/test/unit/solver/test_slv1t.py index 3cbe419..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 @@ -87,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): @@ -100,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): @@ -113,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 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