diff --git a/apps/tutorials/codegen/HeatEquationKernelWithMaterial.py b/apps/tutorials/codegen/HeatEquationKernelWithMaterial.py
index 6348e6356c6452ee94067047be9eea2f8298721a..ea80613df51bf5d4ea953162b37f1711c12b35c8 100644
--- a/apps/tutorials/codegen/HeatEquationKernelWithMaterial.py
+++ b/apps/tutorials/codegen/HeatEquationKernelWithMaterial.py
@@ -1,12 +1,12 @@
import os
import sys
-script_dir = os.path.dirname(__file__)
-walberla_dir = os.path.join(script_dir, '..', '..', '..', '..', 'walberla')
-sys.path.append(walberla_dir)
import sympy as sp
import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
-from src.materials.alloys import Ti6Al4V
+script_dir = os.path.dirname(__file__)
+walberla_dir = os.path.join(script_dir, '..', '..', '..', '..', 'walberla')
+sys.path.append(walberla_dir)
+# from src.materials.alloys import Ti6Al4V
from src.materials.alloys.SS316L import SS316L
from src.materials.assignment_converter import assignment_converter
@@ -27,7 +27,13 @@ with CodeGeneration() as ctx:
Ti6Al4V = SS316L.create_SS316L(u.center())
# Convert assignments to pystencils format
+ print("Print statements")
+ print(Ti6Al4V.thermal_diffusivity)
+ print(Ti6Al4V.thermal_diffusivity.expr)
+ print(Ti6Al4V.thermal_diffusivity.assignments)
subexp, subs = assignment_converter(Ti6Al4V.thermal_diffusivity.assignments)
+ print(subexp)
+ print(subs)
subexp.append(ps.Assignment(thermal_diffusivity, Ti6Al4V.thermal_diffusivity.expr))
diff --git a/src/materials/alloys/Alloy.py b/src/materials/alloys/Alloy.py
index 8b638106d9ec9b090857d2b6a0513d250ad89cc1..43f5fe9edf2e77b4fd3adf8da9c08e02939eb088 100644
--- a/src/materials/alloys/Alloy.py
+++ b/src/materials/alloys/Alloy.py
@@ -22,15 +22,15 @@ class Alloy:
Optional properties:
density (PropertyTypes): Density of the alloy.
+ dynamic_viscosity (PropertyTypes): Dynamic viscosity of the alloy.
heat_capacity (PropertyTypes): Heat capacity of the alloy.
heat_conductivity (PropertyTypes): Heat conductivity of the alloy.
+ kinematic_viscosity (PropertyTypes): Kinematic viscosity of the alloy.
latent_heat_of_fusion (PropertyTypes): Latent heat of fusion of the alloy.
latent_heat_of_vaporization (PropertyTypes): Latent heat of vaporization of the alloy.
+ surface_tension (PropertyTypes): Surface tension of the alloy.
thermal_diffusivity (PropertyTypes): Thermal diffusivity of the alloy.
thermal_expansion_coefficient (PropertyTypes): Thermal expansion coefficient of the alloy.
- dynamic_viscosity (PropertyTypes): Dynamic viscosity of the alloy.
- kinematic_viscosity (PropertyTypes): Kinematic viscosity of the alloy.
- surface_tension (PropertyTypes): Surface tension of the alloy.
"""
elements: List[Element]
composition: ArrayTypes # List of fractions summing to 1.0
@@ -44,15 +44,16 @@ class Alloy:
# Optional properties with default values
density: PropertyTypes = None
+ dynamic_viscosity: PropertyTypes = None
heat_capacity: PropertyTypes = None
heat_conductivity: PropertyTypes = None
+ kinematic_viscosity: PropertyTypes = None
latent_heat_of_fusion: PropertyTypes = None
latent_heat_of_vaporization: PropertyTypes = None
+ surface_tension: PropertyTypes = None
thermal_diffusivity: PropertyTypes = None
thermal_expansion_coefficient: PropertyTypes = None
- dynamic_viscosity: PropertyTypes = None
- kinematic_viscosity: PropertyTypes = None
- surface_tension: PropertyTypes = None
+
def solidification_interval(self) -> Tuple[float, float]:
"""
@@ -65,14 +66,17 @@ class Alloy:
def __post_init__(self) -> None:
"""
- Initializes properties based on elemental composition and validates the composition.
+ Initializes properties based on elemental composition and validates the composition and phase transition temperatures.
Raises:
- ValueError: If the sum of the composition array does not equal 1.
+ ValueError: If the sum of the composition array does not equal 1 or if the solidus temperature is greater than the liquidus temperature.
"""
if not np.isclose(sum(self.composition), 1.0, atol=1e-12):
raise ValueError(f"The sum of the composition array must be 1.0, got {sum(self.composition)}")
+ if self.temperature_solidus > self.temperature_liquidus:
+ raise ValueError("The solidus temperature must be less than or equal to the liquidus temperature.")
+
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)
@@ -109,7 +113,7 @@ if __name__ == '__main__':
# Invalid Property Assignment
try:
- Ti64.heat_conductivity = "invalid_value"
+ Ti64.heat_conductivity = "invalid_value" # type: ignore
except TypeError as e:
print(f"Invalid Property Assignment Test Passed: {e}")
diff --git a/src/materials/alloys/SS316L/SS316L.py b/src/materials/alloys/SS316L/SS316L.py
index c79a425595e1781b202dff9599d4c0b6afbb4e6a..2ca4198d617609156af31a97669753b2dbb2335a 100644
--- a/src/materials/alloys/SS316L/SS316L.py
+++ b/src/materials/alloys/SS316L/SS316L.py
@@ -1,11 +1,13 @@
+import numpy as np
import sympy as sp
from pathlib import Path
from typing import Union
from src.materials.alloys.Alloy import Alloy
from src.materials.elements import Fe, Cr, Mn, Ni
-from src.materials.models import thermal_diffusivity_by_heat_conductivity
+from src.materials.models import thermal_diffusivity_by_heat_conductivity, density_by_thermal_expansion, MaterialPropertyWrapper
from src.materials.data_handler.data_handler import read_data, celsius_to_kelvin, thousand_times
from src.materials.interpolators import interpolate_property
+from src.materials.constants import Constants
def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
@@ -16,7 +18,7 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
T (Union[float, sp.Symbol]): Temperature as a symbolic variable or numeric value.
Returns:
- Tuple[Alloy, Tuple, Tuple, Tuple, Tuple]: Initialized SS316L alloy with physical properties and temperature arrays.
+ Alloy: Initialized SS316L alloy with physical properties.
Notes:
- **Data Units**: All input data should be in SI units. If data is not in SI units, it must be converted. For example:
@@ -32,6 +34,8 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
Example:
If you have temperature data in Celsius and property data in non-SI units, convert the temperature to Kelvin and property values to SI units before using them.
"""
+ if isinstance(T, float) and (T < 0 or T > 3000): # Example valid range: 0 to 3000 K
+ raise ValueError("Invalid temperature. Temperature must be within the valid range (0 to 3000 K).")
# Define the alloy with specific elemental composition and phase transition temperatures
SS316L = Alloy(
elements=[Fe, Cr, Mn, Ni],
@@ -47,7 +51,7 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
# Paths to data files using relative paths
density_data_file_path = str(base_dir / 'density_temperature.txt')
heat_capacity_data_file_path = str(base_dir / 'heat_capacity_temperature.txt')
- heat_conductivity_data_file_path = str(base_dir / 'heat_conductivity_temperature.txt')
+ heat_conductivity_data_file_path = str(base_dir / '..' / 'SS316L' / 'heat_conductivity_temperature.txt')
# Read temperature and material property data from the files
density_temp_array, density_array = read_data(density_data_file_path)
@@ -69,24 +73,78 @@ def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
heat_conductivity_temp_array = celsius_to_kelvin(heat_conductivity_temp_array)
# Perform interpolation for each property
- SS316L.density = interpolate_property(T, density_temp_array, density_array)
+ print('density')
+ # SS316L.density = interpolate_property(T, density_temp_array, density_array)
+ SS316L.density = MaterialPropertyWrapper(density_by_thermal_expansion(
+ T, Constants.temperature_room, 7940.7476, SS316L.thermal_expansion_coefficient))
+ print('SS316L.density: ', SS316L.density)
+ print('heat_capacity')
SS316L.heat_capacity = interpolate_property(T, heat_capacity_temp_array, heat_capacity_array)
+ print('heat_conductivity')
SS316L.heat_conductivity = interpolate_property(T, heat_conductivity_temp_array, heat_conductivity_array)
# Calculate thermal diffusivity from thermal conductivity, density, and heat capacity
- SS316L.thermal_diffusivity = interpolate_property(T, density_temp_array,
+ print('thermal_diffusivity')
+
+ if isinstance(T, sp.Symbol):
+ thermal_diffusivity_array = thermal_diffusivity_by_heat_conductivity(
+ SS316L.heat_conductivity.evalf(T, density_temp_array),
+ SS316L.density.evalf(T, heat_conductivity_temp_array),
+ # SS316L.density.expr,
+ SS316L.heat_capacity.evalf(T, heat_conductivity_temp_array)
+ )
+ print('SS316L.heat_conductivity.evalf(T, density_temp_array): ', (SS316L.heat_conductivity.evalf(T, density_temp_array))) # <class 'numpy.ndarray'>
+ print('SS316L.density.expr: ', type(SS316L.density.expr)) # <class 'sympy.core.mul.Mul'>
+ print('SS316L.heat_capacity.evalf(T, heat_conductivity_temp_array): ', (SS316L.heat_capacity.evalf(T, heat_conductivity_temp_array))) # <class 'numpy.ndarray'>
+ print('thermal_diffusivity_array: ', type(thermal_diffusivity_array)) # <class 'sympy.tensor.array.dense_ndim_array.ImmutableDenseNDimArray'>
+
+ # Calculate thermal diffusivity values for each temperature in the array
+ thermal_diffusivity_array1 = np.array([
+ (thermal_diffusivity_by_heat_conductivity(
+ (SS316L.heat_conductivity.evalf(T, temp)),
+ # SS316L.density.evalf(T, temp),
+ SS316L.density.expr,
+ (SS316L.heat_capacity.evalf(T, temp))
+ ))
+ for temp in density_temp_array
+ ])
+
+ print('thermal_diffusivity_array1 type:', type(thermal_diffusivity_array1))
+ print('thermal_diffusivity_array1 dtype:', thermal_diffusivity_array1.dtype)
+ print('thermal_diffusivity_array1 shape:', thermal_diffusivity_array1.shape)
+ print('First few values of thermal_diffusivity_array1:', thermal_diffusivity_array1[:5])
+
+ SS316L.thermal_diffusivity = interpolate_property(T, density_temp_array, thermal_diffusivity_array)
+
+ elif isinstance(T, float):
+ # elif isinstance(T, sp.Symbol):
+ # Calculate thermal diffusivity from thermal conductivity, density, and heat capacity
+ thermal_diffusivity_array2 = []
+ for temp in heat_conductivity_temp_array:
+ k = SS316L.heat_conductivity.evalf(T, temp)
+ rho = SS316L.density.evalf(T, temp)
+ cp = SS316L.heat_capacity.evalf(T, temp)
+ thermal_diffusivity2 = thermal_diffusivity_by_heat_conductivity(k, rho, cp)
+ thermal_diffusivity_array2.append(thermal_diffusivity2)
+
+ # Interpolate thermal diffusivity using the temperature array and calculated values
+ SS316L.thermal_diffusivity = interpolate_property(T, density_temp_array, thermal_diffusivity_array2)
+
+
+ '''SS316L.thermal_diffusivity = interpolate_property(T, density_temp_array,
thermal_diffusivity_by_heat_conductivity(
SS316L.heat_conductivity.evalf(T, density_temp_array),
- SS316L.density.evalf(T, density_temp_array),
- SS316L.heat_capacity.evalf(T, density_temp_array)
- ))
+ # SS316L.density.evalf(T, heat_conductivity_temp_array),
+ SS316L.density.expr,
+ SS316L.heat_capacity.evalf(T, heat_conductivity_temp_array)
+ ))'''
return SS316L
if __name__ == '__main__':
Temp = sp.Symbol('Temp')
- alloy = create_SS316L(Temp)
+ alloy = create_SS316L(300.15)
# Print the composition of each element in the alloy
for i in range(len(alloy.composition)):
diff --git a/src/materials/alloys/Ti6Al4V.py b/src/materials/alloys/Ti6Al4V.py
index 6ba95809333b53cd0da0f361f9de0438678f8218..49b6d8ec09b9fbd1fd649b7fdaf056c0ee9e3e22 100644
--- a/src/materials/alloys/Ti6Al4V.py
+++ b/src/materials/alloys/Ti6Al4V.py
@@ -4,7 +4,7 @@ from src.materials.alloys.Alloy import Alloy
from src.materials.constants import Constants
from src.materials.elements import Ti, Al, V
from src.materials.interpolators import interpolate_equidistant, interpolate_lookup, interpolate_property
-from src.materials.models import density_by_thermal_expansion, thermal_diffusivity_by_heat_conductivity
+from src.materials.models import density_by_thermal_expansion, thermal_diffusivity_by_heat_conductivity, MaterialPropertyWrapper
from src.materials.typedefs import MaterialProperty
@@ -54,7 +54,7 @@ def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
)
# Compute density based on thermal expansion and other constants
- Ti6Al4V.density = MaterialProperty(density_by_thermal_expansion(
+ Ti6Al4V.density = MaterialPropertyWrapper(density_by_thermal_expansion(
T, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient))
# Interpolate heat conductivity based on temperature
@@ -66,9 +66,16 @@ def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
thermal_diffusivity_by_heat_conductivity(
Ti6Al4V.heat_conductivity.evalf(T, Ti6Al4V.solidification_interval()),
Ti6Al4V.density.expr,
+ # Ti6Al4V.density.evalf(T, Ti6Al4V.solidification_interval()),
Ti6Al4V.heat_capacity
))
+ print('thermal_diffusivity_by_heat_conductivity: ', thermal_diffusivity_by_heat_conductivity(
+ Ti6Al4V.heat_conductivity.evalf(T, Ti6Al4V.solidification_interval()),
+ Ti6Al4V.density.expr,
+ # Ti6Al4V.density.evalf(T, Ti6Al4V.solidification_interval()),
+ Ti6Al4V.heat_capacity
+ ))
return Ti6Al4V
diff --git a/src/materials/assignment_converter.py b/src/materials/assignment_converter.py
index 043634c55363e556ef89c06feb88073a97b52bd1..e397304cad5ad714bbc7799fbae475013e7223fb 100644
--- a/src/materials/assignment_converter.py
+++ b/src/materials/assignment_converter.py
@@ -4,28 +4,43 @@ import pystencils as ps
from typing import List, Dict, Tuple, Union
from pystencils.types.quick import Arr, Fp
from pystencils.sympyextensions.typed_sympy import CastFunc
-from src.materials.typedefs import Assignment
+from src.materials.typedefs import Assignment, ArrayTypes
-def type_mapping(type_str: str) -> Union[np.dtype, Arr]:
+def type_mapping(type_str: str, length: int) -> Union[np.dtype, Arr]:
"""
Maps a string representation of a type to a corresponding numpy or pystencils data type.
Parameters:
type_str (str): The string representation of the type (e.g., "double[]", "float[]").
+ length (int): The length of the array for array types.
Returns:
Union[np.dtype, Arr]: The corresponding numpy or pystencils data type.
- - "double[]" maps to Arr(Fp(64)) which represents a 64-bit floating point array.
- - "float[]" maps to Arr(Fp(32)) which represents a 32-bit floating point array.
- - Other type strings are mapped to corresponding numpy data types.
+
+ Examples:
+ >>> type_mapping("double[]", 5)
+ PsArrayType(element_type=PsIeeeFloatType( width=64, const=True ), size=5, const=False)
+ >>> type_mapping("float[]", 3)
+ PsArrayType(element_type=PsIeeeFloatType( width=32, const=True ), size=3, const=False)
+ >>> type_mapping("double", 1)
+ dtype('float64')
"""
if type_str == "double[]":
- return Arr(Fp(64)) # 64-bit floating point array
+ return Arr(Fp(64, const=True), length=length) # 64-bit floating point array
elif type_str == "float[]":
- return Arr(Fp(32)) # 32-bit floating point array
+ return Arr(Fp(32, const=True), length=length) # 32-bit floating point array
+ elif type_str == "double":
+ return np.dtype('float64')
+ elif type_str == "float":
+ return np.dtype('float32')
+ elif type_str == "int":
+ return np.dtype('int32')
+ elif type_str == "bool":
+ return np.dtype('bool')
else:
- return np.dtype(type_str) # Default to numpy data type
+ # return np.dtype(type_str) # Default to numpy data type
+ raise ValueError(f"Unsupported type string: {type_str}")
def assignment_converter(assignment_list: List[Assignment]) \
@@ -50,7 +65,15 @@ def assignment_converter(assignment_list: List[Assignment]) \
subs = {}
for assignment in assignment_list:
- ps_type = type_mapping(assignment.lhs_type)
+ if assignment.lhs is None or assignment.rhs is None or assignment.lhs_type is None:
+ raise ValueError("Invalid assignment: lhs, rhs, and lhs_type must not be None")
+
+ for assignment in assignment_list:
+ if isinstance(assignment.rhs, ArrayTypes):
+ length = len(assignment.rhs)
+ else:
+ length = 1
+ ps_type = type_mapping(assignment.lhs_type, length)
data_symbol = ps.TypedSymbol(assignment.lhs.name, ps_type)
if isinstance(ps_type, Arr):
@@ -78,9 +101,9 @@ if __name__ == '__main__':
# Test type_mapping function
print("Testing type_mapping function:")
- print(type_mapping("double[]")) # Expected: Arr(Fp(64))
- print(type_mapping("float[]")) # Expected: Arr(Fp(32))
- print(type_mapping("int")) # Expected: np.int32 or equivalent numpy type
+ print(type_mapping("double[]", 5)) # Expected: Arr(Fp(64), length=5)
+ print(type_mapping("float[]", 3)) # Expected: Arr(Fp(32), length=3)
+ print(type_mapping("int", 1)) # Expected: np.int32 or equivalent numpy type
# Test assignment_converter function
print("\nTesting assignment_converter function:")
diff --git a/src/materials/interpolators.py b/src/materials/interpolators.py
index 21b97788136949eccdd92a6b6b66440ed84d8072..9cb496430790a7582c34b62d05077504da682cbd 100644
--- a/src/materials/interpolators.py
+++ b/src/materials/interpolators.py
@@ -1,6 +1,7 @@
import numpy as np
import sympy as sp
-from typing import Union
+from typing import Union, List, Tuple
+from src.materials.models import Wrapper, MaterialPropertyWrapper
from src.materials.typedefs import Assignment, ArrayTypes, MaterialProperty
COUNT = 0
@@ -10,7 +11,8 @@ def interpolate_equidistant(
T: Union[float, sp.Symbol],
T_base: float,
T_incr: float,
- v_array: ArrayTypes) -> Union[float, MaterialProperty]:
+ # v_array: ArrayTypes) -> Union[float, MaterialProperty]: # PropertyTypes
+ v_array: ArrayTypes) -> MaterialProperty:
"""
Perform equidistant interpolation for symbolic or numeric temperature values.
@@ -35,20 +37,20 @@ def interpolate_equidistant(
sym_idx = sp.symbols(label + '_idx', cls=sp.Idx)
sym_dec = sp.Symbol(label + "_dec")
- T_base = sp.Float(T_base)
- T_incr = sp.Float(T_incr)
+ T_base = Wrapper(T_base) # T_base = sp.Float(T_base)
+ T_incr = Wrapper(T_incr) # T_incr = sp.Float(T_incr)
pos = (T - T_base) / T_incr
pos_dec = pos - sp.floor(pos)
result = sp.Piecewise(
- (sp.Float(float(v_array[0])), T < T_base),
- (sp.Float(float(v_array[-1])), T >= T_base + (len(v_array) - 1) * T_incr),
+ (Wrapper(v_array[0]), T < T_base),
+ (Wrapper(v_array[-1]), T >= T_base + (len(v_array) - 1) * T_incr),
((1 - sym_dec) * sym_data[sym_idx] + sym_dec * sym_data[sym_idx + 1], True)
)
-
+ print('sym_data, type(sym_data): ', sym_data, type(sym_data))
sub_expressions = [
- Assignment(sym_data, tuple(sp.Float(float(v)) for v in v_array), "double[]"),
+ Assignment(sym_data, tuple(Wrapper(v) for v in v_array), "double[]"),
Assignment(sym_idx, pos, "int"),
Assignment(sym_dec, pos_dec, "double")
]
@@ -61,16 +63,22 @@ def interpolate_equidistant(
min_temp = T_base
max_temp = T_base + (n - 1) * T_incr
if T < min_temp:
- return float(v_array[0])
+ # return float(v_array[0])
+ # return MaterialProperty(sp.Float(float(v_array[0])), [])
+ return MaterialPropertyWrapper(v_array[0])
elif T > max_temp:
- return float(v_array[-1])
+ # return float(v_array[-1])
+ # return MaterialProperty(sp.Float(float(v_array[-1])), [])
+ return MaterialPropertyWrapper(v_array[-1])
pos = (T - T_base) / T_incr
pos_int = int(pos)
pos_dec = pos - pos_int
value = (1 - pos_dec) * v_array[pos_int] + pos_dec * v_array[pos_int + 1]
- return float(value)
+ # return float(value)
+ # return MaterialProperty(sp.Float(float(value)), [])
+ return MaterialPropertyWrapper(value)
else:
raise ValueError(f"Unsupported type for T: {type(T)}")
@@ -79,7 +87,8 @@ def interpolate_equidistant(
def interpolate_lookup(
T: Union[float, sp.Symbol],
T_array: ArrayTypes,
- v_array: ArrayTypes) -> Union[float, MaterialProperty]:
+ # v_array: ArrayTypes) -> Union[float, MaterialProperty]: # PropertyTypes
+ v_array: ArrayTypes) -> MaterialProperty:
"""
Perform lookup-based interpolation for symbolic or numeric temperature values.
@@ -97,8 +106,8 @@ def interpolate_lookup(
v_array = np.flip(v_array)
if isinstance(T, sp.Symbol):
- T_array = [sp.Float(float(t)) for t in T_array]
- v_array = [v if isinstance(v, sp.Expr) else sp.Float(float(v)) for v in v_array]
+ T_array = [Wrapper(t) for t in T_array]
+ v_array = [Wrapper(v) for v in v_array]
conditions = [
(v_array[0], T < T_array[0]),
@@ -114,12 +123,18 @@ def interpolate_lookup(
T_array = np.asarray(T_array)
v_array = np.asarray(v_array)
if T < T_array[0]:
- return float(v_array[0])
+ # return float(v_array[0])
+ #nreturn MaterialProperty(sp.Float(float(v_array[0])), [])
+ return MaterialPropertyWrapper(v_array[0])
elif T > T_array[-1]:
- return float(v_array[-1])
+ # return float(v_array[-1])
+ # return MaterialProperty(sp.Float(float(v_array[-1])), [])
+ return MaterialPropertyWrapper(v_array[-1])
else:
- return float(np.interp(T, T_array, v_array))
-
+ # return float(np.interp(T, T_array, v_array))
+ # return MaterialProperty(sp.Float(float(np.interp(T, T_array, v_array))), [])
+ v = np.interp(T, T_array, v_array)
+ return MaterialPropertyWrapper(v)
else:
raise ValueError(f"Invalid input for T: {type(T)}")
@@ -141,8 +156,9 @@ def interpolate_property(
T: Union[float, sp.Symbol],
temp_array: ArrayTypes,
prop_array: ArrayTypes,
- temp_array_limit=6,
- force_lookup=False) -> Union[float, MaterialProperty]:
+ temp_array_limit: int = 6,
+ # force_lookup=False) -> Union[float, MaterialProperty]: # PropertyTypes
+ force_lookup: bool = False) -> MaterialProperty:
"""
Perform interpolation based on the type of temperature array (equidistant or not).
@@ -154,6 +170,13 @@ def interpolate_property(
:return: Interpolated value.
:raises ValueError: If T_array and v_array lengths do not match.
"""
+ # Ensure that temp_array and prop_array are arrays
+ if not isinstance(temp_array, (Tuple, List, np.ndarray)):
+ raise TypeError(f"Expected temp_array to be a list or array, got {type(temp_array)}")
+
+ if not isinstance(prop_array, (Tuple, List, np.ndarray)):
+ raise TypeError(f"Expected prop_array to be a list or array, got {type(prop_array)}")
+
if len(temp_array) != len(prop_array):
raise ValueError("T_array and v_array must have the same length")
@@ -164,10 +187,12 @@ def interpolate_property(
incr = test_equidistant(np.asarray(temp_array))
if force_lookup or incr == 0 or len(temp_array) < temp_array_limit:
+ print('interpolate_lookup')
return interpolate_lookup(T, temp_array, prop_array)
else:
+ print('interpolate_equidistant')
return interpolate_equidistant(
- T, float(temp_array[0]), float(incr), prop_array)
+ T, float(temp_array[0]), incr, prop_array)
if __name__ == '__main__':
@@ -183,59 +208,36 @@ if __name__ == '__main__':
# Interpolation Lookup Tests
print("Interpolate Lookup Tests:")
symbolic_result = interpolate_lookup(Temp, density_temp_array1, density_v_array1)
- print("Symbolic interpolation result with interpolate_lookup (arrays):",
- symbolic_result.expr if isinstance(symbolic_result, MaterialProperty) else symbolic_result)
-
+ print("Symbolic interpolation result with interpolate_lookup (arrays):", symbolic_result.expr)
symbolic_result_list = interpolate_lookup(Temp, density_temp_array2, density_v_array2)
- print("Symbolic interpolation result with interpolate_lookup (lists):",
- symbolic_result_list.expr if isinstance(symbolic_result_list, MaterialProperty) else symbolic_result_list)
-
+ print("Symbolic interpolation result with interpolate_lookup (lists):", symbolic_result_list.expr)
numeric_result = interpolate_lookup(1894.7, density_temp_array2, density_v_array2)
- print("Numeric interpolation result with interpolate_lookup (arrays):", numeric_result)
-
- 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)
-
- try:
- error_result = interpolate_lookup("what's up", density_temp_array2, density_v_array2)
- except ValueError as e:
- print("Error:", e)
+ print("Numeric interpolation result with interpolate_lookup (arrays):", numeric_result.expr)
+ 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.expr)
# Interpolation Equidistant Tests
print("\nInterpolate Equidistant Tests:")
symbolic_result_eq = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array1)
print("Symbolic interpolation result with interpolate_equidistant (numpy arrays):", symbolic_result_eq.expr)
- print("Assignments:", symbolic_result_eq.assignments)
-
symbolic_result_eq_list = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array2)
print("Symbolic interpolation result with interpolate_equidistant (lists):", symbolic_result_eq_list.expr)
- print("Assignments:", symbolic_result_eq_list.assignments)
-
numeric_result_eq = interpolate_equidistant(1900.2, Temp_base, Temp_incr, density_v_array1)
- print("Numeric interpolation result with interpolate_equidistant (numpy arrays):", numeric_result_eq)
-
+ print("Numeric interpolation result with interpolate_equidistant (numpy arrays):", numeric_result_eq.expr)
numeric_result_eq_list = interpolate_equidistant(1900.2, Temp_base, Temp_incr, density_v_array2)
- print("Numeric interpolation result with interpolate_equidistant (lists):", numeric_result_eq_list)
-
- try:
- error_result_eq = interpolate_equidistant("what's up", Temp_base, Temp_incr, density_v_array2)
- except ValueError as e:
- print("Error:", e)
+ print("Numeric interpolation result with interpolate_equidistant (lists):", numeric_result_eq_list.expr)
# Consistency Test
print("\nConsistency Test between interpolate_lookup and interpolate_equidistant:")
test_temps = [150.0, 200.0, 350.0, 450.0, 550.0, 650.0]
for TT in test_temps:
- lookup_result = interpolate_lookup(TT, np.flip(np.array([600.0, 500.0, 400.0, 300.0, 200.0])),
- np.array([10.0, 20.0, 30.0, 40.0, 50.0]))
+ lookup_result = interpolate_lookup(TT, np.flip(np.array([600.0, 500.0, 400.0, 300.0, 200.0])), np.array([10.0, 20.0, 30.0, 40.0, 50.0]))
equidistant_result = interpolate_equidistant(TT, 200.0, 100.0, np.array([10.0, 20.0, 30.0, 40.0, 50.0]))
print(f"Temperature {TT}:")
- print(f"interpolate_lookup result: {lookup_result}")
- print(f"interpolate_equidistant result: {equidistant_result}")
-
- if isinstance(lookup_result, (int, float)) and isinstance(equidistant_result, (int, float)):
- assert np.isclose(lookup_result, equidistant_result), f"Inconsistent results for T = {TT}"
+ print(f"interpolate_lookup result: {lookup_result.expr}")
+ print(f"interpolate_equidistant result: {equidistant_result.expr}")
+ if isinstance(lookup_result.expr, (int, float)) and isinstance(equidistant_result.expr, (int, float)):
+ assert np.isclose(lookup_result.expr, equidistant_result.expr), f"Inconsistent results for T = {TT}"
else:
print(f"Skipping comparison for symbolic or non-numeric result at T = {TT}")
diff --git a/src/materials/models.py b/src/materials/models.py
index cc1f1f088b2e9438f5ce29e39d7c5a9151661393..56014a9752a79bf0c0a7c34d461c75808d2fb88d 100644
--- a/src/materials/models.py
+++ b/src/materials/models.py
@@ -1,6 +1,31 @@
import numpy as np
import sympy as sp
-from typing import Union
+from typing import Union, List
+from src.materials.typedefs import MaterialProperty
+
+
+def Wrapper(value: Union[sp.Expr, int, float, np.ndarray]) -> Union[sp.Expr, List[sp.Expr]]:
+ if isinstance(value, sp.Expr):
+ return sp.sympify(value)
+ if isinstance(value, (int, float)):
+ return sp.Float(float(value))
+ if isinstance(value, np.ndarray):
+ return [sp.Float(float(v)) for v in value]
+ # return [Wrapper(v) for v in value]
+ raise ValueError("Unsupported type for value in Wrapper")
+
+
+'''def MaterialPropertyWrapper(value: Union[sp.Expr, int, float, np.ndarray]) -> MaterialProperty:
+ wrapped_value = Wrapper(value)
+ if isinstance(wrapped_value, list):
+ # Instead of returning a list of MaterialProperties, wrap the entire array
+ return MaterialProperty(np.array(wrapped_value))
+ return MaterialProperty(wrapped_value)'''
+
+
+def MaterialPropertyWrapper(value: Union[sp.Expr, int, float, np.ndarray]) -> MaterialProperty:
+ wrapped_value = Wrapper(value)
+ return MaterialProperty(wrapped_value)
def density_by_thermal_expansion(
@@ -24,9 +49,17 @@ def density_by_thermal_expansion(
:param thermal_expansion_coefficient: Thermal expansion coefficient. Can be a numeric value (float), numpy array (np.ndarray), or a symbolic expression (sp.Expr).
:return: Density at temperature T. Returns a float if all inputs are numeric; otherwise, returns a numpy array or a sympy expression.
"""
+ # Input validation
+ for param, name in zip([temperature, temperature_base, density_base, thermal_expansion_coefficient], ['temperature', 'temperature_base', 'density_base', 'thermal_expansion_coefficient']):
+ if isinstance(param, (float, int)) and param <= 0:
+ raise ValueError(f"{name.capitalize()} must be positive.")
+ if isinstance(param, np.ndarray) and (param <= 0).any():
+ raise ValueError(f"{name.capitalize()} array contains non-positive values.")
+
if isinstance(temperature, sp.Symbol) and isinstance(thermal_expansion_coefficient, np.ndarray):
raise ValueError("This combination is not valid. Provide the temperature data which is used for the thermal_expansion_coefficient and create a new interpolation material property.")
- return density_base * (1 + thermal_expansion_coefficient * (temperature - temperature_base)) ** (-3)
+ result = density_base * (1 + thermal_expansion_coefficient * (temperature - temperature_base)) ** (-3)
+ return result
def thermal_diffusivity_by_heat_conductivity(
@@ -49,8 +82,16 @@ def thermal_diffusivity_by_heat_conductivity(
:param heat_capacity: Specific heat capacity. Can be a numeric value (float), numpy array (np.ndarray), or a symbolic expression (sp.Expr).
:return: Thermal diffusivity. Returns a float if all inputs are numeric; otherwise, returns a numpy array or a sympy expression.
"""
+ # Input validation
+ for param, name in zip([heat_conductivity, density, heat_capacity], ['heat_conductivity', 'density', 'heat_capacity']):
+ if isinstance(param, (float, int)) and param <= 0:
+ raise ValueError(f"{name.capitalize()} must be positive.")
+ if isinstance(param, np.ndarray) and (param <= 0).any():
+ raise ValueError(f"{name.capitalize()} array contains non-positive values.")
+
if isinstance(density, sp.Piecewise) or isinstance(heat_conductivity, sp.Piecewise) or isinstance(heat_capacity, sp.Piecewise):
raise ValueError("Do not insert piecewise functions here. Prepare the data with corresponding temperature intervals.")
+
result = heat_conductivity / (density * heat_capacity)
return result
@@ -100,4 +141,86 @@ if __name__ == '__main__':
c_p_symbolic = sp.Symbol('c_p')
rho_symbolic = sp.Symbol('rho')
print(f"Test 6 - Symbolic thermal diffusivity computation: {thermal_diffusivity_by_heat_conductivity(k_symbolic, rho_symbolic, c_p_symbolic)}")
- # Expected output: Sympy expression for thermal diffusivity as a function of k, rho, and c_p
+
+ # New test case for mixed input types
+ print("\nTest 9: Mixed input types")
+
+ # For density_by_thermal_expansion
+ T_mixed = np.linspace(800, 1000, 3) # numpy array
+ T_base = 293.15 # float
+ rho_base = 8000.0 # float
+ tec_symbolic = sp.Symbol('alpha') # symbolic
+
+ try:
+ result_density = density_by_thermal_expansion(T_mixed, T_base, rho_base, tec_symbolic)
+ print(f"Density with mixed inputs: {result_density}")
+ except Exception as e:
+ print(f"Error in density calculation with mixed inputs: {str(e)}")
+
+ # For thermal_diffusivity_by_heat_conductivity
+ k_mixed = np.linspace(20, 30, 3) # numpy array
+ rho_symbolic = sp.Symbol('rho') # symbolic
+ cp_float = 500.0 # float
+
+ try:
+ result_diffusivity = thermal_diffusivity_by_heat_conductivity(k_mixed, rho_symbolic, cp_float)
+ print(f"Thermal diffusivity with mixed inputs: {result_diffusivity}")
+ except Exception as e:
+ print(f"Error in thermal diffusivity calculation with mixed inputs: {str(e)}")
+
+ # Test case for all symbolic inputs
+ print("\nTest 10: All symbolic inputs")
+ T_sym = sp.Symbol('T')
+ tec_sym = sp.Symbol('alpha')
+ k_sym = sp.Symbol('k')
+ rho_sym = sp.Symbol('rho')
+ cp_sym = sp.Symbol('cp')
+
+ density_sym = density_by_thermal_expansion(T_sym, T_base, rho_base, tec_sym)
+ print(f"Symbolic density: {density_sym}")
+
+ diffusivity_sym = thermal_diffusivity_by_heat_conductivity(k_sym, rho_sym, cp_sym)
+ print(f"Symbolic thermal diffusivity: {diffusivity_sym}")
+
+ # Test 11: Edge cases with zero and negative values
+ print("\nTest 11: Edge cases with zero and negative values")
+ try:
+ print(density_by_thermal_expansion(-273.15, T_0, rho, tec)) # Absolute zero
+ except ValueError as e:
+ print(f"Handled error for negative temperature: {e}")
+
+ try:
+ print(density_by_thermal_expansion(T_1, T_0, rho, -1e-6)) # Negative TEC
+ except ValueError as e:
+ print(f"Handled error for negative TEC: {e}")
+
+ # Test 12: Large arrays
+ print("\nTest 12: Large arrays")
+ large_T = np.linspace(1000, 2000, 10000)
+ large_k = np.linspace(30, 40, 10000)
+ try:
+ large_density = density_by_thermal_expansion(large_T, T_0, rho, tec)
+ print(f"Large density array: {large_density}")
+ except Exception as e:
+ print(f"Error with large density array: {e}")
+
+ try:
+ large_diffusivity = thermal_diffusivity_by_heat_conductivity(large_k, rho, c_p)
+ print(f"Large diffusivity array: {large_diffusivity}")
+ except Exception as e:
+ print(f"Error with large diffusivity array: {e}")
+
+ # Test 13: Complex numbers (if applicable)
+ print("\nTest 13: Complex numbers")
+ complex_T = sp.Symbol('T') + sp.I
+ try:
+ complex_density = density_by_thermal_expansion(complex_T, T_0, rho, tec)
+ print(f"Complex density: {complex_density}")
+ except Exception as e:
+ print(f"Error with complex density: {e}")
+
+ try:
+ complex_diffusivity = thermal_diffusivity_by_heat_conductivity(k_symbolic, rho_symbolic + sp.I, c_p_symbolic)
+ print(f"Complex diffusivity: {complex_diffusivity}")
+ except Exception as e:
+ print(f"Error with complex diffusivity: {e}")
diff --git a/src/materials/tests.py b/src/materials/tests.py
new file mode 100644
index 0000000000000000000000000000000000000000..40bdbc08f555522de9412e9def731d9ed6692f08
--- /dev/null
+++ b/src/materials/tests.py
@@ -0,0 +1,326 @@
+import pytest
+import numpy as np
+import sympy as sp
+from pystencils.types.quick import Arr
+
+from src.materials.alloys.Alloy import Alloy
+from src.materials.alloys.SS316L.SS316L import create_SS316L
+from src.materials.assignment_converter import type_mapping, assignment_converter
+from src.materials.elements import Fe, Cr, Mn, Ni
+from src.materials.interpolators import interpolate_property, interpolate_lookup, interpolate_equidistant, \
+ test_equidistant
+from src.materials.models import density_by_thermal_expansion, thermal_diffusivity_by_heat_conductivity
+from src.materials.typedefs import MaterialProperty, Assignment
+
+
+def test_alloy_creation():
+ # Test creating an alloy with valid elemental composition and phase transition temperatures
+ alloy = Alloy(elements=[Fe, Cr, Mn, Ni], composition=[0.7, 0.2, 0.05, 0.05], temperature_solidus=1700, temperature_liquidus=1800)
+ assert np.allclose(alloy.composition, [0.7, 0.2, 0.05, 0.05])
+ assert alloy.temperature_solidus == 1700
+ assert alloy.temperature_liquidus == 1800
+
+ # Test creating an alloy with invalid elemental composition
+ with pytest.raises(ValueError):
+ Alloy(elements=[Fe, Cr], composition=[0.6, 0.5], temperature_solidus=1700, temperature_liquidus=1800)
+
+ # Test creating an alloy with invalid phase transition temperatures
+ with pytest.raises(ValueError):
+ Alloy(elements=[Fe, Cr, Mn, Ni], composition=[0.7, 0.2, 0.05, 0.05], temperature_solidus=1900, temperature_liquidus=1800)
+
+def test_create_SS316L():
+ # Test creating SS316L alloy with a float temperature input
+ alloy = create_SS316L(1400.0)
+ assert isinstance(alloy.density, MaterialProperty)
+ assert isinstance(alloy.heat_capacity, MaterialProperty)
+ assert isinstance(alloy.heat_conductivity, MaterialProperty)
+ assert isinstance(alloy.thermal_diffusivity, MaterialProperty)
+
+ # Test creating SS316L alloy with a symbolic temperature input
+ T = sp.Symbol('T')
+ alloy = create_SS316L(T)
+ assert isinstance(alloy.density, MaterialProperty)
+ assert isinstance(alloy.heat_capacity, MaterialProperty)
+ assert isinstance(alloy.heat_conductivity, MaterialProperty)
+ assert isinstance(alloy.thermal_diffusivity, MaterialProperty)
+
+ # Test creating SS316L alloy with different combinations of property calculation methods
+ alloy = create_SS316L(1400.0)
+ alloy.density = alloy.density.expr
+ alloy.heat_capacity = alloy.heat_capacity.evalf(T, 1400.0)
+ assert isinstance(alloy.density, sp.Expr)
+ assert isinstance(alloy.heat_capacity, float)
+
+def test_interpolate_property():
+ # Test interpolating properties with float temperature inputs
+ T_symbol = sp.Symbol('T') # Define a symbolic variable
+ T_value = 1400.0
+ temp_array = np.array([1300.0, 1400.0, 1500.0])
+ prop_array = np.array([100.0, 200.0, 300.0])
+ result = interpolate_property(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test interpolating properties with symbolic temperature inputs
+ result = interpolate_property(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test interpolating properties with different combinations of temperature and property arrays
+ temp_array = [1300.0, 1400.0, 1500.0]
+ prop_array = [100.0, 200.0, 300.0]
+ result = interpolate_property(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ temp_array = tuple([1300.0, 1400.0, 1500.0])
+ prop_array = tuple([100.0, 200.0, 300.0])
+ result = interpolate_property(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test with ascending and descending order arrays
+ temp_array_asc = np.array([1300.0, 1400.0, 1500.0])
+ prop_array_asc = np.array([100.0, 200.0, 300.0])
+ result_asc = interpolate_property(T_symbol, temp_array_asc, prop_array_asc)
+ assert isinstance(result_asc, MaterialProperty)
+ assert np.isclose(result_asc.evalf(T_symbol, T_value), 200.0)
+
+ temp_array_desc = np.array([1500.0, 1400.0, 1300.0])
+ prop_array_desc = np.array([300.0, 200.0, 100.0])
+ result_desc = interpolate_property(T_symbol, temp_array_desc, prop_array_desc)
+ assert isinstance(result_desc, MaterialProperty)
+ assert np.isclose(result_desc.evalf(T_symbol, T_value), 200.0)
+
+ # Test with different temp_array_limit values
+ result_limit_6 = interpolate_property(T_symbol, temp_array_asc, prop_array_asc, temp_array_limit=6)
+ assert isinstance(result_limit_6, MaterialProperty)
+ assert np.isclose(result_limit_6.evalf(T_symbol, T_value), 200.0)
+
+ result_limit_3 = interpolate_property(T_symbol, temp_array_asc, prop_array_asc, temp_array_limit=3)
+ assert isinstance(result_limit_3, MaterialProperty)
+ assert np.isclose(result_limit_3.evalf(T_symbol, T_value), 200.0)
+
+ # Test with force_lookup True and False
+ result_force_lookup_true = interpolate_property(T_symbol, temp_array_asc, prop_array_asc, force_lookup=True)
+ assert isinstance(result_force_lookup_true, MaterialProperty)
+ assert np.isclose(result_force_lookup_true.evalf(T_symbol, T_value), 200.0)
+
+ result_force_lookup_false = interpolate_property(T_symbol, temp_array_asc, prop_array_asc, force_lookup=False)
+ assert isinstance(result_force_lookup_false, MaterialProperty)
+ assert np.isclose(result_force_lookup_false.evalf(T_symbol, T_value), 200.0)
+
+def test_interpolate_lookup():
+ # Define a symbolic variable
+ T_symbol = sp.Symbol('T')
+
+ # Test lookup interpolation with float temperature inputs
+ T_value = 1400.0
+ temp_array = np.array([1300.0, 1400.0, 1500.0])
+ prop_array = np.array([100.0, 200.0, 300.0])
+ result = interpolate_lookup(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test lookup interpolation with symbolic temperature inputs
+ result = interpolate_lookup(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test with different combinations of temperature and property arrays
+ temp_array = [1300.0, 1400.0, 1500.0]
+ prop_array = [100.0, 200.0, 300.0]
+ result = interpolate_lookup(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ temp_array = tuple([1300.0, 1400.0, 1500.0])
+ prop_array = tuple([100.0, 200.0, 300.0])
+ result = interpolate_lookup(T_symbol, temp_array, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test with ascending and descending order arrays
+ temp_array_asc = np.array([1300.0, 1400.0, 1500.0])
+ prop_array_asc = np.array([100.0, 200.0, 300.0])
+ result_asc = interpolate_lookup(T_symbol, temp_array_asc, prop_array_asc)
+ assert isinstance(result_asc, MaterialProperty)
+ assert np.isclose(result_asc.evalf(T_symbol, T_value), 200.0)
+
+ temp_array_desc = np.array([1500.0, 1400.0, 1300.0])
+ prop_array_desc = np.array([300.0, 200.0, 100.0])
+ result_desc = interpolate_lookup(T_symbol, temp_array_desc, prop_array_desc)
+ assert isinstance(result_desc, MaterialProperty)
+ assert np.isclose(result_desc.evalf(T_symbol, T_value), 200.0)
+
+def test_interpolate_equidistant():
+ # Define a symbolic variable
+ T_symbol = sp.Symbol('T')
+
+ # Test equidistant interpolation with float temperature inputs
+ T_value = 1400.0
+ temp_base = 1300.0
+ temp_incr = 100.0
+ prop_array = np.array([100.0, 200.0, 300.0])
+ result = interpolate_equidistant(T_value, temp_base, temp_incr, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test equidistant interpolation with symbolic temperature inputs
+ result = interpolate_equidistant(T_symbol, temp_base, temp_incr, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test with different combinations of temperature and property arrays
+ temp_base = 1300.0
+ temp_incr = 100.0
+ prop_array = [100.0, 200.0, 300.0]
+ result = interpolate_equidistant(T_symbol, temp_base, temp_incr, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ temp_base = 1300.0
+ temp_incr = 100.0
+ prop_array = tuple([100.0, 200.0, 300.0])
+ result = interpolate_equidistant(T_symbol, temp_base, temp_incr, prop_array)
+ assert isinstance(result, MaterialProperty)
+ assert np.isclose(result.evalf(T_symbol, T_value), 200.0)
+
+ # Test with ascending and descending order arrays
+ temp_base = 1300.0
+ temp_incr = 100.0
+ prop_array_asc = np.array([100.0, 200.0, 300.0])
+ result_asc = interpolate_equidistant(T_symbol, temp_base, temp_incr, prop_array_asc)
+ assert isinstance(result_asc, MaterialProperty)
+ assert np.isclose(result_asc.evalf(T_symbol, T_value), 200.0)
+
+ temp_base = 1300.0
+ temp_incr = 100.0
+ prop_array_desc = np.array([300.0, 200.0, 100.0])
+ result_desc = interpolate_equidistant(T_symbol, temp_base, temp_incr, prop_array_desc)
+ assert isinstance(result_desc, MaterialProperty)
+ assert np.isclose(result_desc.evalf(T_symbol, T_value), 200.0)
+
+# Define the temp fixture
+@pytest.fixture
+def temp():
+ return np.array([100.0, 200.0, 300.0, 400.0, 500.0])
+
+def test_test_equidistant(temp):
+ # Test with equidistant temperature array
+ assert test_equidistant(temp) == 100.0
+
+ # Test with non-equidistant temperature array
+ temp_non_equidistant = np.array([100.0, 150.0, 300.0, 450.0, 500.0])
+ assert test_equidistant(temp_non_equidistant) == 0.0
+
+def test_density_by_thermal_expansion():
+ # Test calculating density with float temperature and thermal expansion coefficient inputs
+ T = 1400.0
+ T_ref = 1000.0
+ rho_ref = 8000.0
+ alpha = 1e-5
+ result = density_by_thermal_expansion(T, T_ref, rho_ref, alpha)
+ assert np.isclose(result, 7904.76)
+
+ # Test calculating density with symbolic temperature and thermal expansion coefficient inputs
+ T = sp.Symbol('T')
+ result = density_by_thermal_expansion(T, T_ref, rho_ref, alpha)
+ assert isinstance(result, sp.Expr)
+ assert np.isclose(result.subs(T, 1400.0), 7904.76)
+
+ # Test calculating density with numpy array temperature and thermal expansion coefficient inputs
+ T = np.array([1300.0, 1400.0, 1500.0])
+ result = density_by_thermal_expansion(T, T_ref, rho_ref, alpha)
+ assert isinstance(result, np.ndarray)
+ assert np.allclose(result, [7928.43, 7904.76, 7881.19])
+
+def test_thermal_diffusivity_by_heat_conductivity():
+ # Test calculating thermal diffusivity with float heat conductivity, density, and heat capacity inputs
+ k = 50.0
+ rho = 8000.0
+ c_p = 500.0
+ result = thermal_diffusivity_by_heat_conductivity(k, rho, c_p)
+ assert np.isclose(result, 1.25e-5)
+
+ # Test calculating thermal diffusivity with symbolic heat conductivity, density, and heat capacity inputs
+ k = sp.Symbol('k')
+ rho = sp.Symbol('rho')
+ c_p = sp.Symbol('c_p')
+ result = thermal_diffusivity_by_heat_conductivity(k, rho, c_p)
+ assert isinstance(result, sp.Expr)
+ assert result == k / (rho * c_p)
+
+ # Test calculating thermal diffusivity with numpy array heat conductivity, density, and heat capacity inputs
+ k = np.array([50.0, 60.0, 70.0])
+ rho = np.array([8000.0, 8100.0, 8200.0])
+ c_p = np.array([500.0, 510.0, 520.0])
+ result = thermal_diffusivity_by_heat_conductivity(k, rho, c_p)
+ assert isinstance(result, np.ndarray)
+ assert np.allclose(result, [1.25e-5, 1.45243282e-5, 1.64165103e-5])
+
+def test_alloy_single_element():
+ # Test creating an alloy with a single element
+ alloy = Alloy(elements=[Fe], composition=[1.0], temperature_solidus=1800, temperature_liquidus=1900)
+ assert len(alloy.elements) == 1
+ assert alloy.composition == [1.0]
+
+def test_alloy_property_modification():
+ # Test accessing and modifying individual alloy properties
+ alloy = create_SS316L(1400.0)
+ assert isinstance(alloy.density, MaterialProperty)
+ # Set the density to a MaterialProperty instance with a constant expression
+ alloy.density = MaterialProperty(expr=sp.Float(8000.0))
+ assert alloy.density.expr == sp.Float(8000.0)
+
+def test_create_SS316L_invalid_temperature():
+ # Test creating SS316L alloy with invalid temperature inputs
+ with pytest.raises(ValueError):
+ create_SS316L(-100.0) # Negative temperature
+ with pytest.raises(ValueError):
+ create_SS316L(3500.0) # Temperature outside valid range
+
+def test_density_by_thermal_expansion_invalid_inputs():
+ # Test calculating density with invalid temperature or thermal expansion coefficient inputs
+ with pytest.raises(ValueError):
+ density_by_thermal_expansion(-100.0, 1000.0, 8000.0, 1e-5)
+ with pytest.raises(ValueError):
+ density_by_thermal_expansion(1000.0, -1000.0, 8000.0, 1e-5)
+ with pytest.raises(ValueError):
+ density_by_thermal_expansion(1000.0, 1000.0, -8000.0, 1e-5)
+ with pytest.raises(ValueError):
+ density_by_thermal_expansion(1000.0, 1000.0, 8000.0, -1e-5)
+ with pytest.raises(ValueError):
+ density_by_thermal_expansion(np.array([-100.0, 1000.0]), 1000.0, 8000.0, 1e-5)
+
+def test_thermal_diffusivity_invalid_inputs():
+ # Test calculating thermal diffusivity with invalid heat conductivity, density, or heat capacity inputs
+ with pytest.raises(ValueError):
+ thermal_diffusivity_by_heat_conductivity(-10.0, 8000.0, 500.0) # Negative heat conductivity
+ with pytest.raises(ValueError):
+ thermal_diffusivity_by_heat_conductivity(50.0, -8000.0, 500.0) # Negative density
+ with pytest.raises(ValueError):
+ thermal_diffusivity_by_heat_conductivity(50.0, 8000.0, -500.0) # Negative heat capacity
+
+def test_type_mapping_unsupported():
+ # Test mapping unsupported type strings
+ with pytest.raises(ValueError):
+ type_mapping("unsupported_type", 1)
+
+# Additional test cases
+def test_type_mapping():
+ assert isinstance(type_mapping("double[]", 5), Arr)
+ assert isinstance(type_mapping("float[]", 3), Arr)
+ assert type_mapping("double", 1) == np.dtype('float64')
+ assert type_mapping("float", 1) == np.dtype('float32')
+ assert type_mapping("int", 1) == np.dtype('int32')
+ assert type_mapping("bool", 1) == np.dtype('bool')
+
+ with pytest.raises(ValueError):
+ type_mapping("unsupported_type", 1)
+
+def test_assignment_converter_invalid():
+ # Test converting assignments with invalid or missing attributes
+ invalid_assignment = Assignment(lhs=None, rhs=None, lhs_type=None)
+ with pytest.raises(ValueError, match="Invalid assignment: lhs, rhs, and lhs_type must not be None"):
+ assignment_converter([invalid_assignment])
diff --git a/src/materials/typedefs.py b/src/materials/typedefs.py
index 7ba5ff70446c44c87c2cbf31ca62c22eae29ea84..a51c4f26202fec27775c1e3a93621f335c59f5ab 100644
--- a/src/materials/typedefs.py
+++ b/src/materials/typedefs.py
@@ -43,7 +43,7 @@ class Assignment:
# Type aliases for various data types used in the project
ArrayTypes = Union[np.ndarray, List, Tuple] # Numerical values can be represented as numpy arrays, lists, or tuples
-
+# def evalf(self, symbol: sp.Symbol, temperature: Union[sp.Symbol, float, ArrayTypes]) -> Union[float, np.ndarray]:
@dataclass
class MaterialProperty:
@@ -78,10 +78,8 @@ class MaterialProperty:
if not self.expr.free_symbols:
return float(self.expr)
- # If temperature is a numpy array, list, or tuple, evaluate the property for each temperature
+ # If temperature is a numpy array, list, or tuple (ArrayTypes), evaluate the property for each temperature
if isinstance(temperature, get_args(ArrayTypes)):
- if isinstance(temperature, np.ndarray):
- temperature = temperature.astype(float) # Ensure all values are Python floats
return np.array([self.evalf(symbol, t) for t in temperature])
# Convert any numpy scalar to Python float
@@ -101,8 +99,12 @@ class MaterialProperty:
# Evaluate the material property with the substitutions
result = sp.N(self.expr.subs(substitutions))
- return float(result)
-
+ # return float(result)
+ # Try to convert the result to float if possible
+ try:
+ return float(result)
+ except TypeError:
+ return result # Return the symbolic expression if it cannot be converted to a float
# Type alias for properties that can be either constant floats or MaterialProperty instances
PropertyTypes = Union[float, MaterialProperty]
@@ -195,3 +197,59 @@ if __name__ == '__main__':
# Test 10: Property with no symbolic dependencies
mp7 = MaterialProperty(sp.Float(1500.))
print(f"Test 10 - Property with no symbolic dependencies: {mp7.evalf(T, 200.)}") # Expected output: 1500.0
+
+ # Test 11: Piecewise function
+ mp8 = MaterialProperty(sp.Piecewise((100, T < 0), (200, T >= 100), (T, True)))
+ print(f"Test 11a - Piecewise function at T=-10: {mp8.evalf(T, -10)}") # Expected: 100
+ print(f"Test 11b - Piecewise function at T=50: {mp8.evalf(T, 50)}") # Expected: 50
+ print(f"Test 11c - Piecewise function at T=150: {mp8.evalf(T, 150)}") # Expected: 200
+
+ # Test 12: Complex expression with trigonometric functions
+ mp9 = MaterialProperty(100 * sp.sin(T) + 50 * sp.cos(T))
+ print(f"Test 12a - Complex expression at T=0: {mp9.evalf(T, 0)}") # Expected: 50
+ print(f"Test 12b - Complex expression at T=pi/2: {mp9.evalf(T, np.pi/2)}") # Expected: 100
+
+ # Test 13: Expression with logarithmic and exponential functions
+ mp10 = MaterialProperty(10 * sp.log(T + 1) + 5 * sp.exp(T/100))
+ print(f"Test 13 - Log and exp expression at T=10: {mp10.evalf(T, 10)}") # Expected: ~28.19
+
+ # Test 14: Property with multiple variables
+ T2 = sp.Symbol('T2')
+ mp11 = MaterialProperty(T * T2)
+ res = mp11.evalf(T, 10)
+ print(f"Test 14 - Multi-variable property: {res}") # This should raise an error or return a symbolic expression
+ assert isinstance(res, sp.Expr)
+ assert sp.Eq(res, 10 * T2)
+
+ # Test 15: Handling very large and very small numbers
+ mp12 = MaterialProperty(1e20 * T + 1e-20)
+ print(f"Test 15a - Large numbers: {mp12.evalf(T, 1e5)}")
+ print(f"Test 15b - Small numbers: {mp12.evalf(T, 1e-5)}")
+
+ # Test 16: Property with rational functions
+ mp13 = MaterialProperty((T**2 + 1) / (T + 2))
+ print(f"Test 16 - Rational function at T=3: {mp13.evalf(T, 3)}")
+
+ # Test 17: Handling undefined values (division by zero)
+ mp14 = MaterialProperty(1 / (T - 1))
+ try:
+ res = mp14.evalf(T, 1)
+ print(f"Test 17 - Division by zero at T=1: {res}") # "zoo" in SymPy represents complex infinity
+ except ZeroDivisionError:
+ print("Test 17 - Division by zero at T=1: ZeroDivisionError raised as expected")
+
+ # Test 18: Property with absolute value
+ mp15 = MaterialProperty(sp.Abs(T - 50))
+ print(f"Test 18a - Absolute value at T=30: {mp15.evalf(T, 30)}")
+ print(f"Test 18b - Absolute value at T=70: {mp15.evalf(T, 70)}")
+
+ # Test 19: Property with floor and ceiling functions
+ mp16 = MaterialProperty(sp.floor(T/10) + sp.ceiling(T/10))
+ print(f"Test 19 - Floor and ceiling at T=25: {mp16.evalf(T, 25)}")
+
+ # Test 20: Handling complex numbers
+ mp17 = MaterialProperty(sp.sqrt(T))
+ print(f"Test 20a - Square root at T=4: {mp17.evalf(T, 4)}")
+ res = mp17.evalf(T, -1)
+ print(f"Test 20b - Square root at T=-1: {res}")
+ assert isinstance(res, sp.Expr)