diff --git a/apps/tutorials/codegen/HeatEquationKernel.py b/apps/tutorials/codegen/HeatEquationKernel.py
index 8339fdff1cca5859ca71b8e8d0e5c06601c7c1f2..1ce85747a63f7a0dd3146ea11f22030a884a3cbe 100644
--- a/apps/tutorials/codegen/HeatEquationKernel.py
+++ b/apps/tutorials/codegen/HeatEquationKernel.py
@@ -1,6 +1,7 @@
import sys
sys.path.append('/local/ca00xebo/repos/walberla')
import sympy as sp
+import numpy as np
import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
from src.materials.alloys import Ti6Al4V
@@ -20,15 +21,47 @@ with CodeGeneration() as ctx:
heat_pde_discretized = discretize(heat_pde)
heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
Ti6Al4V = Ti6Al4V.create_Ti6Al4V(u.center())
- # density, density_subexpr = Ti6Al4V.density # Alternative
subexp = []
# Iterate over the dictionary items
+ # {'density': [(density, array([50., 40., 30., 20., 10.])), (density_idx, floor(Temp/100) - 2)]}
for label, assignment_list in Ti6Al4V.assignments.items():
- for data_label, data_values in assignment_list:
- data_symbol = ps.TypedSymbol(data_label, Arr(Fp(64))) # Define the symbol with appropriate type
- subexp.append(ps.Assignment(data_symbol, data_values))
- # thermal_diffusivity_expr = thermal_diffusivity_by_heat_conductivity(Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
- thermal_diffusivity_expr = thermal_diffusivity_by_heat_conductivity(Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity) # Alternative
+ print('assignment_list: ', assignment_list)
+
+ '''# Convert to a list of SymPy floats
+ values_as_list = [sp.Float(float(v)) for v in assignment_list[0][1]]
+ print('values_as_list: ', values_as_list)
+ print('type(values_as_list): ', type(values_as_list))'''
+
+ '''data_label = assignment_list[0][0]
+ data_values = assignment_list[0][1]'''
+
+ for assignment in assignment_list:
+ '''print('data_label: ', data_label)
+ print('type(data_label): ', type(data_label))
+ print('data_values: ', data_values)
+ print('type(data_values): ', type(data_values))'''
+
+ # values_as_tuples = tuple(sp.Float(float(v)) for v in data_values)
+
+ if isinstance(assignment.rhs, tuple):
+ data_symbol = ps.TypedSymbol(assignment.lhs.name, Arr(Fp(64)))
+ subexp.append(ps.Assignment(data_symbol, assignment.rhs))
+ elif isinstance(assignment.rhs, sp.Expr):
+ data_symbol = ps.TypedSymbol(assignment.lhs.name, dtype='int64')
+ subexp.append(ps.Assignment(data_symbol, assignment.rhs))
+
+ # data_symbol_array = ps.TypedSymbol(assignment_list[0][0].label.name, data_type='list', dtype='float64')
+ # data_symbol_array = ps.TypedSymbol(assignment_list[0][0].name, ps.types.PsArrayType(ps.create_type("float64")))
+ # data_symbol = ps.TypedSymbol(assignment_list[0][0].label.name, Arr(Fp(64))) #
+
+ # data_symbol_array = sp.IndexedBase(data_symbol_array, shape=len(assignment_list[0][1]))
+
+ # subexp.append(ps.Assignment(data_symbol, values_as_list)) #
+
+ # data_symbol_idx = ps.TypedSymbol(assignment_list[1][0].name, dtype='int64')
+ # subexp.append(ps.Assignment(data_symbol_idx, (assignment_list[1][1]))) # maybe the CastFunc to Int is necessary here
+
+ thermal_diffusivity_expr = thermal_diffusivity_by_heat_conductivity(Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
subexp.append(ps.Assignment(thermal_diffusivity, thermal_diffusivity_expr))
ac = ps.AssignmentCollection(
subexpressions=subexp,
@@ -43,4 +76,6 @@ with CodeGeneration() as ctx:
# ac = ps.simp.simplifications.add_subexpressions_for_field_reads(ac)
# ac = ps.simp.sympy_cse(ac)
+ print(ac)
+
generate_sweep(ctx, 'HeatEquationKernel', ac, varying_parameters=((data_type, str(thermal_diffusivity)),))
diff --git a/src/materials/alloys/Alloy.py b/src/materials/alloys/Alloy.py
index c5ff7250bb6d7b830f95957845d951a3483cdf00..00c9e681ea9a15b97f1317700c15210ea0d66d5c 100644
--- a/src/materials/alloys/Alloy.py
+++ b/src/materials/alloys/Alloy.py
@@ -1,16 +1,15 @@
import numpy as np
-from sympy import Expr
-from typing import List, Union, Dict, Tuple, Any
-from dataclasses import dataclass, field, asdict
-from src.materials.elements import Element
-from src.materials.elements import Ti, Al, V
+from dataclasses import dataclass, field
+from src.materials.elements import Element, Ti, Al, V
+from src.materials.elements import interpolate_atomic_number, interpolate_atomic_mass, interpolate_temperature_boil
+from src.materials.typedefs import ValueTypes, AssignmentType, PropertyTypes
# Alloy Base Class
@dataclass
class Alloy:
- elements: List[Element]
- composition: List[float]
+ elements: list[Element]
+ composition: ValueTypes
atomic_number: float = field(init=False)
atomic_mass: float = field(init=False)
@@ -19,34 +18,33 @@ class Alloy:
temperature_solidus: float
temperature_liquidus: float
- latent_heat_of_fusion: Union[float, Expr] = None
- latent_heat_of_vaporization: Union[float, Expr] = None
- density: Union[float, Expr] = None
- thermal_expansion_coefficient: Union[float, Expr] = None
- heat_capacity: Union[float, Expr] = None
- heat_conductivity: Union[float, Expr] = None
- thermal_diffusivity: Union[float, Expr] = None
+ latent_heat_of_fusion: PropertyTypes = None
+ latent_heat_of_vaporization: PropertyTypes = None
+ density: PropertyTypes = None
+ thermal_expansion_coefficient: PropertyTypes = None
+ heat_capacity: PropertyTypes = None
+ heat_conductivity: PropertyTypes = None
+ thermal_diffusivity: PropertyTypes = None
- # assignments: list[tuple, ...] = field(default_factory=list, init=False)
- assignments: Dict[str, List[Tuple[str, Any]]] = field(default_factory=dict, init=False)
+ assignments: AssignmentType = field(default_factory=dict, init=False)
- def interpolate(self, quantity: str) -> float:
- comp = np.asarray(self.composition)
- qnt = np.asarray([element.__getattribute__(quantity) for element in self.elements])
- return float(np.sum(comp * qnt))
+ def solidification_interval(self) -> [float, float]:
+ return [self.temperature_solidus, self.temperature_liquidus]
def __post_init__(self):
- self.atomic_number = self.interpolate("atomic_number")
- self.atomic_mass = self.interpolate("atomic_mass")
- self.temperature_boil = self.interpolate("temperature_boil")
+ # if sum(self.composition) != 1.:
+ # raise ValueError("The sum of the composition array has to be 1 got", sum(self.composition))
+ if not np.isclose(sum(self.composition), 1.0, atol=1e-8):
+ raise ValueError("The sum of the composition array must be 1, got", sum(self.composition))
- def solidification_interval(self) -> [float, float]:
- return [self.temperature_solidus, self.temperature_liquidus]
+ self.atomic_number = interpolate_atomic_number(self.elements, self.composition)
+ self.atomic_mass = interpolate_atomic_mass(self.elements, self.composition)
+ self.temperature_boil = interpolate_temperature_boil(self.elements, self.composition)
if __name__ == '__main__':
Ti64 = Alloy([Ti, Al, V], [0.90, 0.06, 0.04], 1878, 1928)
- print(0.9*22+0.06*13+0.04*23, Ti64.interpolate("atomic_number"), Ti64.atomic_number)
+ print(0.9 * 22 + 0.06 * 13 + 0.04 * 23, Ti64.atomic_number)
print(Ti64.heat_conductivity)
Ti64.heat_conductivity = 34
print(Ti64.heat_conductivity)
diff --git a/src/materials/alloys/Ti6Al4V.py b/src/materials/alloys/Ti6Al4V.py
index 52a0026ad197f03fc2cab26f42370b639f833fc6..56e2fa0d19ec224a95185a9c04bbe3dd31ad9bdd 100644
--- a/src/materials/alloys/Ti6Al4V.py
+++ b/src/materials/alloys/Ti6Al4V.py
@@ -1,14 +1,13 @@
import sympy as sp
-from typing import Union
from src.materials.alloys.Alloy import Alloy
from src.materials.elements import Ti, Al, V
from src.materials.constants import Constants
from src.materials import interpolators
from src.materials import models
-# from src.materials.property import Property
+from src.materials.typedefs import VariableTypes, AssignmentType
-def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
+def create_Ti6Al4V(T: VariableTypes) -> Alloy:
# 90% Ti, 6% Al, 4% V
Ti6Al4V = Alloy(
[Ti, Al, V],
@@ -18,41 +17,62 @@ def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
)
Ti6Al4V.latent_heat_of_fusion = 290000.0 # m²/s²
- Ti6Al4V.density = 4430.0 # kg/m³ Ref [5]
Ti6Al4V.thermal_expansion_coefficient = 10e-6 # /K Ref[6]
- heat_capacity_expr, heat_capacity_assignments = interpolators.interpolate_equidistant(
- T, 1820, 20, [700, 800, 900, 1000, 1100, 1200], 'heat_capacity')
- Ti6Al4V.heat_capacity = heat_capacity_expr
+ # ---------- density ----------
- # Merge the returned assignments into the Ti6Al4V assignments dictionary
- for label, data in heat_capacity_assignments.items():
- if label in Ti6Al4V.assignments:
- Ti6Al4V.assignments[label].extend(data)
- else:
- Ti6Al4V.assignments[label] = data
-
- Ti6Al4V.heat_conductivity = interpolators.interpolate_lookup(
- T, Ti6Al4V.solidification_interval(), [21.58, 33.60]) # W/m-k
+ # Ti6Al4V.density = 4430.0 # kg/m³ Ref [5]
'''Ti6Al4V.density = models.density_by_thermal_expansion(
- T, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient)
+ T, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient)'''
+
+ # Using interpolate_equidistant for density
+ density_result, sub_expr_density = interpolators.interpolate_equidistant(
+ T, 1878, 8, [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], 'density')
+ Ti6Al4V.density = density_result
+ Ti6Al4V.assignments.update({"density": sub_expr_density})
+
+ # older code kept for reference (look at interpolate_equidistant for this implementation)
+ '''Ti6Al4V.density = interpolators.interpolate_equidistant(
+ T, 1878, 8, [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], Ti6Al4V.assignments, 'density')'''
+ # Error: pystencils.backend.kernelcreation.freeze.FreezeError: Don't know how to freeze expression density_idx
+
+ # Using interpolate_lookup for heat conductivity
+ '''Ti6Al4V.density = interpolators.interpolate_lookup(
+ T, Ti6Al4V.solidification_interval(), [4236.3, 4227.])'''
+
+ # Ti6Al4V.density = interpolators.interpolate_lookup(
+ # T, [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0],
+ # [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
+
+ # ---------- heat_capacity ----------
+
+ # Ti6Al4V.heat_capacity = 526.
- Ti6Al4V.density = interpolators.interpolate_lookup(
- T, [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0],
- [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
+ # Using interpolate_equidistant for heat capacity
+ '''heat_capacity_result, sub_expr_heat_capacity = interpolators.interpolate_equidistant(
+ T, 1820, 20, [526, 540, 550, 560, 570, 580], 'heat_capacity')
+ Ti6Al4V.heat_capacity = heat_capacity_result
+ Ti6Al4V.assignments.update({"heat_capacity": sub_expr_heat_capacity})'''
- density_expr, density_assignments = interpolators.interpolate_equidistant(
- T, 1878.0, 5,
- [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], 'density')
- Ti6Al4V.density = density_expr
+ # Using interpolate_lookup for heat conductivity
+ Ti6Al4V.heat_capacity = interpolators.interpolate_lookup(
+ T, Ti6Al4V.solidification_interval(), [526, 580.0])
- # Merge the returned assignments into the Ti6Al4V assignments dictionary
- for label, data in density_assignments.items():
- if label in Ti6Al4V.assignments:
- Ti6Al4V.assignments[label].extend(data)
- else:
- Ti6Al4V.assignments[label] = data'''
+ # ---------- heat_conductivity ----------
+
+ # Ti6Al4V.heat_conductivity = 6.7
+
+ # Using interpolate_equidistant for heat conductivity
+ '''Ti6Al4V.heat_conductivity, sub_expr_heat_conductivity = interpolators.interpolate_equidistant(
+ T, 20, 1, [6.7, 6.5, 6.3, 6.1, 5.9, 5.7], 'heat_conductivity')
+ Ti6Al4V.assignments.update({"heat_conductivity": sub_expr_heat_conductivity})'''
+
+ # Using interpolate_lookup for heat conductivity
+ Ti6Al4V.heat_conductivity = interpolators.interpolate_lookup(
+ T, Ti6Al4V.solidification_interval(), [21.58, 33.60])
+
+ # ---------- thermal_diffusivity ----------
Ti6Al4V.thermal_diffusivity = models.thermal_diffusivity_by_heat_conductivity(
Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
@@ -63,10 +83,45 @@ def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
if __name__ == '__main__':
Temp = sp.Symbol('Temp')
alloy = create_Ti6Al4V(Temp)
- print('Ti6Al4V')
- for field in vars(alloy):
- print(field, "=", alloy.__getattribute__(field))
- print(models.thermal_diffusivity_by_heat_conductivity(500, 5000, 600))
-
- print(alloy.density.subs(Temp, 1890))
+ print("Testing Ti6Al4V with symbolic temperature:")
+ print('Ti6Al4V Properties:')
+ for field in vars(alloy):
+ print(f"{field} = {alloy.__getattribute__(field)}")
+
+ # Test interpolate_equidistant
+ print("\nTesting interpolate_equidistant:")
+ test_temp = 1850.0 # Example temperature value
+ result, _ = interpolators.interpolate_equidistant(
+ test_temp, 1820, 20.0, [700, 800, 900, 1000, 1100, 1200], 'heat_capacity')
+ print(f"Interpolated heat capacity at T={test_temp} using interpolate_equidistant: {result}")
+
+ _, assignments = interpolators.interpolate_equidistant(
+ Temp, 1820, 20.0, [700, 800, 900, 1000, 1100, 1200], 'heat_capacity')
+ print(f"Assignments: {assignments}")
+
+ result_density, _ = interpolators.interpolate_equidistant(
+ test_temp, 1878, 8, [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], 'density')
+ print(f"Interpolated density at T={test_temp} using interpolate_equidistant: {result_density}")
+
+ # Test interpolate_lookup
+ print("\nTesting interpolate_lookup:")
+ test_T_array = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+ test_v_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+ result_lookup = interpolators.interpolate_lookup(
+ Temp, test_T_array, test_v_array)
+ print(f"Interpolated value using interpolate_lookup: {result_lookup}")
+
+ # Test thermal diffusivity calculation
+ heat_conductivity = 500 # Example value for testing
+ density = 5000 # Example value for testing
+ heat_capacity = 600 # Example value for testing
+ thermal_diffusivity = models.thermal_diffusivity_by_heat_conductivity(
+ heat_conductivity, density, heat_capacity)
+ print(f"Calculated thermal diffusivity: {thermal_diffusivity}")
+
+ # sp.plot(alloy.heat_conductivity, (Temp, 1800, 2000))
+
+ # sp.plot(alloy.density, (Temp, 1800, 2000))
+
+ # sp.plot(alloy.heat_capacity.subs([(alloy.assignments['heat_capacity'][1][0], alloy.assignments['heat_capacity'][1][1]), (alloy.assignments['heat_capacity'][0][0], alloy.assignments['heat_capacity'][0][1])]), (Temp, 1800, 2000))
diff --git a/src/materials/constants.py b/src/materials/constants.py
index 9c804c57558a058feee8e95d2923c5bad65c8131..a1caa883ca6ccca913d92a1ec261d2ec64d1ed8a 100644
--- a/src/materials/constants.py
+++ b/src/materials/constants.py
@@ -3,6 +3,6 @@ from dataclasses import dataclass
@dataclass
class Constants:
- temperature_room = 298.15 # k
- N_a = 6.022141e23 # /mol
- u = 1.660538e-27 # kg
\ No newline at end of file
+ temperature_room = 298.15 # k
+ N_a = 6.022141e23 # /mol
+ u = 1.660538e-27 # kg
diff --git a/src/materials/elements.py b/src/materials/elements.py
index 049fc1756d98550a2f56042807f88124cc66608b..17cf2e7c19664c374243911df163cf2fdf075c48 100644
--- a/src/materials/elements.py
+++ b/src/materials/elements.py
@@ -1,3 +1,4 @@
+import numpy as np
from dataclasses import dataclass
from src.materials.constants import Constants
@@ -41,3 +42,24 @@ V = Element(
latent_heat_of_fusion=21500,
latent_heat_of_vaporization=444000
)
+
+
+def interpolate(values, composition) -> float:
+ c = np.asarray(composition)
+ v = np.asarray(values)
+ return float(np.sum(c * v))
+
+
+def interpolate_atomic_number(elements, composition) -> float:
+ values = [element.atomic_number for element in elements]
+ return interpolate(values, composition)
+
+
+def interpolate_temperature_boil(elements, composition) -> float:
+ values = [element.temperature_boil for element in elements]
+ return interpolate(values, composition)
+
+
+def interpolate_atomic_mass(elements, composition) -> float:
+ values = [element.atomic_mass for element in elements]
+ return interpolate(values, composition)
diff --git a/src/materials/interpolators.py b/src/materials/interpolators.py
index 8d57a4bf0bb339952dd62c77dcf910c654b0f5f6..8b945f6df614ffc2831e66dd561d0471d3d1e356 100644
--- a/src/materials/interpolators.py
+++ b/src/materials/interpolators.py
@@ -1,127 +1,232 @@
import numpy as np
import sympy as sp
-from typing import Union, List, Tuple, Dict
-# import pystencils as ps
-# from pystencils.types.quick import *
-# from src.materials.property import Property
-from pystencils.sympyextensions import CastFunc
-
-
-def interpolate_equidistant(T: Union[float, sp.Symbol], T_base: float, T_incr: float,
- data_array: Union[np.ndarray, List[float]], label: str) \
- -> Union[float, Tuple[Union[sp.Expr, float], Dict[str, List[Tuple[str, Tuple[float, ...]]]]]]: # TODO
- """
- Interpolates values for equidistant temperature data points.
-
- :param label:
- :param T: Temperature value, either symbolic (sp.Expr) or numeric (float).
- :param T_base: Base temperature value.
- :param T_incr: Increment between successive temperature data points.
- :param data_array: Array of data values corresponding to equidistant temperature points.
- :return: Interpolated value at temperature T.
- """
+from src.materials.typedefs import ValueTypes, PropertyTypes, VariableTypes, Assignment, AssignmentType
+
+
+# older code kept for reference
+# def interpolate_equidistant(T: VariableTypes, T_base: float, T_incr: float, v_array: ValueTypes, assignments: AssignmentType, label: str) -> PropertyTypes:
+def interpolate_equidistant(T: VariableTypes, T_base: float, T_incr: float, v_array: ValueTypes,
+ label: str) -> PropertyTypes:
if isinstance(T, sp.Symbol):
- # Ensure all data_array elements are compatible with SymPy
- # data_array = [sp.Float(float(x)) for x in data_array]
- # data_symbol = sp.IndexedBase('data_symbol', shape=len(data_array))
- # data_symbol = sp.IndexedBase(ps.TypedSymbol('density_data_symbol', Ptr(Fp(32))), shape=len(data_array))
- data_array = tuple(sp.Float(float(x)) for x in data_array) # Alternative
- # data_symbol = ps.TypedSymbol(label, Arr(Fp(64))) # Alternative
- data_base = sp.IndexedBase(label, shape=len(data_array)) # Alternative
- print(type(data_base))
- print(type(data_base.label.name))
- print(type(data_array))
- # idx = sp.Symbol('idx', integer=True)
- idx = sp.floor((T - T_base) / T_incr)
- idx_int = CastFunc(idx, np.int64)
- print(idx.is_Integer)
- mod = (T - T_base) - idx * T_incr # mod replaces (T - T_base) % T_incr
+ data_base = sp.IndexedBase(label, shape=len(v_array))
+ idx = sp.symbols(label + '_idx', cls=sp.Idx)
+
+ pos = (T - T_base) / T_incr
+ pos_int = sp.floor(pos)
+ pos_dec = sp.frac(pos)
+
result = sp.Piecewise(
- (data_array[0], T < T_base),
- (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
- # (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True))
- (((T_incr - mod) * data_base[idx_int] + (mod * data_base[idx_int + 1])) / T_incr, True)) # Alternative
- # return result.subs(data_symbol, sp.Array(data_array))
- assignments = {label: [(label, data_array)]}
- return result, assignments # [result, ps.Assignment(data_symbol, data_array), ps.Assignment(data_symbol, data_array)] # Alternative resurn a property, 1st arg: result, 2nd: List[tuple]
- else:
+ (sp.Float(float(v_array[0])), T < T_base),
+ (sp.Float(float(v_array[-1])), T >= T_base + (len(v_array) - 1) * T_incr),
+ ((1 - pos_dec) * data_base[idx] + pos_dec * data_base[idx + 1], True))
+
+ # older code kept for reference
+ '''if label in assignments:
+ raise ValueError("Assigning a label again", label)
+ assignments[label] = (Assignment(data_base, tuple([float(v) for v in v_array])), Assignment(idx, pos_int), )'''
+
+ sub_expressions = (
+ Assignment(data_base, tuple(sp.Float(float(v)) for v in v_array)),
+ Assignment(idx, pos_int),
+ )
+
+ # return result # older code kept for reference
+ return result, sub_expressions
+ elif isinstance(T, float):
if T < T_base:
- return data_array[0]
- if T > T_base + (len(data_array) - 1) * T_incr:
- return data_array[-1]
- idx = int((T - T_base) / T_incr)
- mod = (T - T_base) - idx * T_incr # mod replaces (T - T_base) % T_incr
- return ((T_incr - mod) * data_array[idx] + (mod * data_array[idx + 1])) / T_incr
-
-
-def interpolate_lookup(T: Union[float, sp.Symbol], base_temperature: Union[np.ndarray, List[float]],
- base_values: Union[np.ndarray, List[float]]) -> Union[float, sp.Expr]: # TODO
- """
- Interpolates properties based on temperature for both symbolic and numeric cases.
-
- :param T: Temperature value, either symbolic or numeric.
- :param base_temperature: List or array of base temperatures.
- :param base_values: List or array of property values corresponding to base temperatures.
- """
- # Ensure that base_temperature and base_values are numpy arrays or lists of floats
- if isinstance(base_temperature, (np.ndarray, list)):
- base_temperature = np.array(base_temperature, dtype=float)
- if isinstance(base_values, (np.ndarray, list)):
- base_values = np.array(base_values, dtype=float)
-
- # Check if base_temperature and base_values have the same length
- if len(base_temperature) != len(base_values):
- raise ValueError("base_temperature and base_values must have the same length")
- # assert len(base_temperature) == len(base_values), "base_temperature and base_values must have the same length"
+ return v_array[0], ()
+
+ pos = (T - T_base) / T_incr
+ pos_int = int(pos)
+ pos_dec = np.modf(pos)[0]
+
+ if pos_int + 1 >= len(v_array):
+ return v_array[-1], ()
+
+ value = (1 - pos_dec) * v_array[pos_int] + pos_dec * v_array[pos_int + 1]
+ return value, ()
+ else:
+ raise ValueError(f"Unsupported type for T: {type(T)}")
+
+
+def interpolate_lookup(T: VariableTypes, T_array: ValueTypes, v_array: ValueTypes) -> PropertyTypes:
+ if len(T_array) != len(v_array): # Check if T_array and v_array have the same length
+ raise ValueError("T_array and v_array must have the same length")
+ # assert len(T_array) == len(v_array), "T_array and v_array must have the same length"
if isinstance(T, sp.Symbol): # Symbolic expression
- base_temperature = [sp.Float(float(t)) for t in base_temperature]
- base_values = [sp.Float(float(val)) for val in base_values]
+ T_array = tuple([sp.Float(float(t)) for t in T_array]) # is it better to have a list or tuple here?
+ print(T_array)
+ v_array = tuple([sp.Float(float(v)) for v in v_array])
+ print(v_array)
+ # T_array = np.asarray(T_array, dtype=sp.Float)
+ # v_array = np.asarray(v_array, dtype=sp.Float)
+
conditions = [
- (base_values[0], T < base_temperature[0]),
- (base_values[-1], T >= base_temperature[-1])]
- for i in range(len(base_temperature) - 1):
- interp_expr = (base_values[i] +
- (base_values[i + 1] - base_values[i]) / (base_temperature[i + 1] - base_temperature[i]) *
- (T - base_temperature[i]))
- if i + 1 == len(base_temperature) - 1:
- conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
- else:
- conditions.append((interp_expr, True))
- # return sp.Piecewise(*conditions)
+ (v_array[0], T < T_array[0]),
+ (v_array[-1], T >= T_array[-1])]
+ for i in range(len(T_array) - 1):
+ interp_expr = (
+ v_array[i] + (v_array[i + 1] - v_array[i]) / (T_array[i + 1] - T_array[i]) * (T - T_array[i]))
+ conditions.append(
+ (interp_expr, sp.And(T >= T_array[i], T < T_array[i + 1]) if i + 1 < len(T_array) - 1 else True))
return sp.Piecewise(*conditions)
- else: # Numeric expression
- try:
- T = float(T) # Ensure T is a float for numeric case
- except (TypeError, ValueError) as e:
- raise ValueError(f"Invalid numeric input for T: {T}") from e
-
- return np.interp(T, base_temperature, base_values)
+ elif isinstance(T, float): # Numeric expression
+ return np.interp(T, T_array, v_array)
+ raise ValueError(f"Invalid input for T: Union[float, sp.Symbol] got", type(T))
+'''
if __name__ == '__main__':
# Example usage:
Temp = sp.Symbol('Temp')
density_temp_array1 = np.array([200.0, 300.0, 400.0, 500.0, 600.0])
- density_data_array1 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
+ density_v_array1 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
density_temp_array2 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
- density_data_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+ density_v_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+
+ assignments = dict()
# Symbolic interpolation
- symbolic_result = interpolate_lookup(Temp, density_temp_array2, density_data_array2)
+ symbolic_result = interpolate_lookup(Temp, density_temp_array2, density_v_array2)
print("Symbolic interpolation result with interpolate_lookup:", symbolic_result)
# Numeric interpolation
- numeric_result = interpolate_lookup(450.0, density_temp_array2, density_data_array2)
+ numeric_result = interpolate_lookup(1900.2, density_temp_array2, density_v_array2)
print("Numeric interpolation result with interpolate_lookup:", numeric_result)
+ # Handle error case for invalid type
+ try:
+ error_result = interpolate_lookup("what's up", density_temp_array2, density_v_array2)
+ except ValueError as e:
+ print("Error:", e)
+
T_base = 200
T_incr = 100
# Symbolic interpolation
- symbolic_result = interpolate_equidistant(Temp, T_base, T_incr, density_data_array2, 'density')
+ symbolic_result, ass1 = interpolate_equidistant(Temp, T_base, T_incr, density_v_array1, 'density')
print("Symbolic interpolation result with interpolate_equidistant:", symbolic_result)
+ print(ass1)
# Numeric interpolation
- numeric_result = interpolate_equidistant(450.0, T_base, T_incr, density_data_array2, 'density')
+ numeric_result = interpolate_equidistant(500.0, T_base, T_incr, density_v_array1, 'density')
print("Numeric interpolation result with interpolate_equidistant:", numeric_result)
+
+ # Ensure proper access to elements in the assignments dictionary
+ for label, assignment_list in assignments.items():
+
+ # Substitute idx and base values into the expression
+ evaluated_result = symbolic_result.subs([
+ (assignments['density'][1][0], assignments['density'][1][1]), # idx substitution
+ (Temp, 450.0), # Temperature substitution
+ (assignments['density'][0][0], assignments['density'][0][1])
+ ])
+ print('evaluated_result: ', evaluated_result)
+
+ # Convert NumPy array to standard Python floats before substitution
+ indexed_data = [sp.Float(float(v)) for v in assignments['density'][0][1]] # Convert to SymPy floats
+ print(indexed_data)
+ indexed_base = assignments['density'][0][0] # The IndexedBase symbol
+ print(indexed_base)
+
+ # map the IndexedBase to the data array for evaluation
+ for i, value in enumerate(indexed_data):
+ evaluated_result = evaluated_result.subs(indexed_base[i], value)
+
+ print(evaluated_result)
+ print(evaluated_result.evalf())
+'''
+
+if __name__ == '__main__':
+ Temp = sp.Symbol('Temp')
+ density_temp_array1 = np.array([200.0, 300.0, 400.0, 500.0, 600.0])
+ density_v_array1 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
+ density_temp_array2 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+ density_v_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+
+ T_base = 200
+ T_incr = 100
+
+ # Example usage of interpolate_lookup
+ print("\nInterpolate Lookup Tests:")
+
+ # Symbolic interpolation with arrays
+ symbolic_result = interpolate_lookup(Temp, density_temp_array1, density_v_array1)
+ print("Symbolic interpolation result with interpolate_lookup (arrays):", symbolic_result)
+
+ # Symbolic interpolation with lists
+ symbolic_result_list = interpolate_lookup(Temp, density_temp_array2, density_v_array2)
+ print("Symbolic interpolation result with interpolate_lookup (lists):", symbolic_result_list)
+
+ # Numeric interpolation with arrays
+ numeric_result = interpolate_lookup(1900.2, density_temp_array2, density_v_array2)
+ print("Numeric interpolation result with interpolate_lookup (arrays):", numeric_result)
+
+ # Numeric interpolation with lists
+ numeric_result_list = interpolate_lookup(1900.2, [200.0, 300.0, 400.0, 500.0, 600.0],
+ [50.0, 40.0, 30.0, 20.0, 10.0])
+ print("Numeric interpolation result with interpolate_lookup (lists):", numeric_result_list)
+
+ # Error handling test for invalid T type
+ try:
+ error_result = interpolate_lookup("what's up", density_temp_array2, density_v_array2)
+ except ValueError as e:
+ print("Error:", e)
+
+ # Example usage of interpolate_equidistant
+ print("\nInterpolate Equidistant Tests:")
+
+ # Symbolic interpolation with numpy arrays
+ symbolic_result_eq, ass1 = interpolate_equidistant(Temp, T_base, T_incr, density_v_array1, 'density')
+ print("Symbolic interpolation result with interpolate_equidistant (numpy arrays):", symbolic_result_eq)
+ print("Assignments:", ass1)
+
+ # Symbolic interpolation with lists
+ symbolic_result_eq_list, ass2 = interpolate_equidistant(Temp, T_base, T_incr, density_v_array2, 'density')
+ print("Symbolic interpolation result with interpolate_equidistant (lists):", symbolic_result_eq_list)
+ print("Assignments:", ass2)
+
+ # Numeric interpolation with numpy arrays
+ numeric_result_eq, _ = interpolate_equidistant(1900.2, T_base, T_incr, density_v_array1, 'density')
+ print("Numeric interpolation result with interpolate_equidistant (numpy arrays):", numeric_result_eq)
+
+ # Numeric interpolation with lists
+ numeric_result_eq_list, _ = interpolate_equidistant(1900.2, T_base, T_incr, density_v_array2, 'density')
+ print("Numeric interpolation result with interpolate_equidistant (lists):", numeric_result_eq_list)
+
+ # Error handling test for invalid T type in interpolate_equidistant
+ try:
+ error_result_eq = interpolate_equidistant("invalid", T_base, T_incr, density_v_array2, 'density')
+ except ValueError as e:
+ print("Error:", e)
+
+ print("\nConsistency Test between interpolate_lookup and interpolate_equidistant:")
+
+ # Consistency check for a range of temperatures
+ test_temps = [150.0, 200.0, 350.0, 450.0, 550.0, 650.0]
+ for T in test_temps:
+ # Numeric results from interpolate_lookup
+ lookup_result = interpolate_lookup(T, density_temp_array1, density_v_array1)
+
+ # Numeric results from interpolate_equidistant
+ equidistant_result, _ = interpolate_equidistant(T, T_base, T_incr, density_v_array1, 'density')
+
+ # Print and assert for consistency
+ print(f"Temperature {T}:")
+ print(f"interpolate_lookup result: {lookup_result}")
+ print(f"interpolate_equidistant result: {equidistant_result}")
+
+ # Allow a small tolerance for floating-point comparison
+ assert np.isclose(lookup_result, equidistant_result), f"Inconsistent results for T = {T}"
+
+ # Consistency check for symbolic interpolation (equivalent results comparison)
+ symbolic_lookup_result = interpolate_lookup(Temp, density_temp_array1, density_v_array1)
+ symbolic_equidistant_result, assignments = interpolate_equidistant(Temp, T_base, T_incr, density_v_array1,
+ 'density')
+
+ # Print symbolic results
+ print("\nSymbolic interpolation results for consistency check:")
+ print(f"interpolate_lookup (symbolic): {symbolic_lookup_result}")
+ print(f"interpolate_equidistant (symbolic): {symbolic_equidistant_result}")
+ print(f"interpolate_equidistant (symbolic): {assignments}")
diff --git a/src/materials/models.py b/src/materials/models.py
index a4a7372eae51d65e9d8454623148c032d962f480..fb6908013a5101ecda23a93cb93fb69cad1b288a 100644
--- a/src/materials/models.py
+++ b/src/materials/models.py
@@ -1,13 +1,13 @@
-from sympy import Expr, Symbol, simplify
-from typing import Union
+import sympy as sp
+from src.materials.typedefs import PropertyTypes, VariableTypes
def density_by_thermal_expansion(
- T: Union[float, Symbol],
+ T: VariableTypes,
base_temperature: float,
base_density: float,
- thermal_expansion_coefficient: Union[float, Expr]) \
- -> Union[float, Expr]:
+ thermal_expansion_coefficient: PropertyTypes) \
+ -> PropertyTypes:
"""
Calculate density based on thermal expansion: rho_0 / (1 + tec * (T - T_0))
@@ -20,10 +20,10 @@ def density_by_thermal_expansion(
def thermal_diffusivity_by_heat_conductivity(
- heat_conductivity: Union[float, Expr],
- density: Union[float, Expr],
- heat_capacity: Union[float, Expr]) \
- -> Union[float, Expr]:
+ heat_conductivity: PropertyTypes,
+ density: PropertyTypes,
+ heat_capacity: PropertyTypes) \
+ -> PropertyTypes:
"""
Calculate thermal diffusivity: k(T) / (rho(T) * c_p(T))
@@ -32,6 +32,6 @@ def thermal_diffusivity_by_heat_conductivity(
:param heat_capacity: Specific heat capacity.
"""
result = heat_conductivity / (density * heat_capacity)
- if isinstance(heat_conductivity, Expr) or isinstance(density, Expr) or isinstance(heat_capacity, Expr):
- result = simplify(result)
+ if isinstance(heat_conductivity, sp.Expr) or isinstance(density, sp.Expr) or isinstance(heat_capacity, sp.Expr):
+ result = sp.simplify(result)
return result
diff --git a/src/materials/property.py b/src/materials/property.py
deleted file mode 100644
index 7a61e3e21e7ed630642153fc77662c631cce3e05..0000000000000000000000000000000000000000
--- a/src/materials/property.py
+++ /dev/null
@@ -1,8 +0,0 @@
-"""from sympy import Expr
-from dataclasses import dataclass
-
-
-@dataclass
-class Property:
- expression: Expr | float
- assignments: tuple[tuple, ...] = tuple"""
diff --git a/src/materials/typedefs.py b/src/materials/typedefs.py
new file mode 100644
index 0000000000000000000000000000000000000000..d56c6d1e7bb6db9b4f1ee95985b7e74018cae474
--- /dev/null
+++ b/src/materials/typedefs.py
@@ -0,0 +1,18 @@
+import numpy as np
+import sympy as sp
+from dataclasses import dataclass
+
+
+@dataclass
+class Assignment:
+ lhs: sp.Symbol
+ rhs: tuple | sp.Expr
+
+
+ValueTypes = np.ndarray | list[float]
+
+PropertyTypes = float | sp.Expr
+
+VariableTypes = float | sp.Symbol
+
+AssignmentType = dict[str, tuple[Assignment, ...]]