diff --git a/apps/tutorials/codegen/HeatEquationKernel.py b/apps/tutorials/codegen/HeatEquationKernel.py
index 6348b0d583e1b5a51052dea7d2c9010731104bae..60a9d5a7e93bff7889bfddfaa6b360a7f2764e9e 100644
--- a/apps/tutorials/codegen/HeatEquationKernel.py
+++ b/apps/tutorials/codegen/HeatEquationKernel.py
@@ -5,7 +5,6 @@ import pystencils as ps
from pystencils_walberla import CodeGeneration, generate_sweep
from src.materials.alloys import Ti6Al4V
from src.materials.alloys.SS316L import SS316L
-from src.materials.models import thermal_diffusivity_by_heat_conductivity
from src.materials.assignment_converter import assignment_converter
with CodeGeneration() as ctx:
@@ -25,13 +24,9 @@ with CodeGeneration() as ctx:
Ti6Al4V = SS316L.create_SS316L(u.center())
# Convert assignments to pystencils format
- subexp, subs = assignment_converter(Ti6Al4V.assignments.items())
+ subexp, subs = assignment_converter(Ti6Al4V.thermal_diffusivity.assignments)
- if Ti6Al4V.thermal_diffusivity is not None:
- subexp.append(ps.Assignment(thermal_diffusivity, Ti6Al4V.thermal_diffusivity))
- else:
- 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))
+ subexp.append(ps.Assignment(thermal_diffusivity, Ti6Al4V.thermal_diffusivity.expr))
ac = ps.AssignmentCollection(
subexpressions=subexp,
@@ -42,6 +37,8 @@ with CodeGeneration() as ctx:
ac = ac.new_with_substitutions(subs)
+
+
# ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
# ac = ps.simp.simplifications.add_subexpressions_for_field_reads(ac)
# ac = ps.simp.sympy_cse(ac)
diff --git a/src/materials/alloys/Alloy.py b/src/materials/alloys/Alloy.py
index cbe3ab2ffbbe62efb7caa6d43ced37da46be9870..23726c14ece0eb4a20bf249637be4b0dff1b19a2 100644
--- a/src/materials/alloys/Alloy.py
+++ b/src/materials/alloys/Alloy.py
@@ -1,18 +1,39 @@
-from typing import List, Union, Tuple
import numpy as np
+from typing import List, Tuple, Union
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, VariableTypes
-from src.materials.interpolators import interpolate_equidistant, interpolate_lookup
-from src.materials.data_handler.data_handler import test_equidistant
+from src.materials.typedefs import ArrayTypes, PropertyTypes
-# Base class for alloys, storing elemental composition and physical properties
@dataclass
class Alloy:
+ """
+ Represents an alloy composed of various elements with given fractions.
+
+ Attributes:
+ elements (List[Element]): List of `Element` objects that make up the alloy.
+ composition (ArrayTypes): List of fractions of each element in the alloy. Fractions should sum to 1.0.
+ temperature_solidus (float): Solidus temperature of the alloy in Kelvin.
+ temperature_liquidus (float): Liquidus temperature of the alloy in Kelvin.
+ atomic_number (float): Average atomic number of the alloy, calculated post-init.
+ atomic_mass (float): Average atomic mass of the alloy, calculated post-init.
+ temperature_boil (float): Boiling temperature of the alloy, calculated post-init.
+
+ Optional properties:
+ density (PropertyTypes): Density of the alloy.
+ heat_capacity (PropertyTypes): Heat capacity of the alloy.
+ heat_conductivity (PropertyTypes): Heat conductivity 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.
+ 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: ValueTypes # List of fractions summing to 1.0
+ composition: ArrayTypes # List of fractions summing to 1.0
temperature_solidus: float
temperature_liquidus: float
@@ -33,65 +54,73 @@ class Alloy:
kinematic_viscosity: PropertyTypes = None
surface_tension: PropertyTypes = None
- # Dictionary for storing calculated assignments
- assignments: AssignmentType = field(default_factory=dict, init=False)
-
def solidification_interval(self) -> Tuple[float, float]:
"""
- Calculate the solidification interval based on solidus and liquidus temperatures.
+ Calculate the solidification interval based on solidus and liquidus temperature of the alloy.
Returns:
- tuple: A tuple containing the solidus and liquidus temperatures.
+ Tuple[float, float]: A tuple containing the solidus and liquidus temperatures.
"""
return self.temperature_solidus, self.temperature_liquidus
- def __post_init__(self):
+ def __post_init__(self) -> None:
"""
- Initialize properties based on elemental composition.
- Ensures the composition sums to 1.0, allowing for floating-point precision issues.
+ Initializes properties based on elemental composition and validates the composition.
+
+ Raises:
+ ValueError: If the sum of the composition array does not equal 1.
"""
if not np.isclose(sum(self.composition), 1.0, atol=1e-12):
- raise ValueError("The sum of the composition array must be 1, got", sum(self.composition))
+ raise ValueError(f"The sum of the composition array must be 1.0, got {sum(self.composition)}")
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)
- def interpolate_property(self, T: VariableTypes, temp_array: ValueTypes, prop_array: ValueTypes, prop_name: str,
- temp_array_limit=6, force_lookup=False):
- """
- Perform interpolation based on the type of temperature array (equidistant or not).
-
- :param T: Temperature, either symbolic (sp.Symbol) or numeric (float).
- :param temp_array: Array of temperatures.
- :param prop_array: Array of property values corresponding to the temperatures.
- :param prop_name: Property name for symbolic representation.
- :param temp_array_limit: Minimum length of the temperature array to force equidistant interpolation.
- :param force_lookup: Boolean to force lookup interpolation regardless of temperature array type.
- :return: Interpolated value.
- :raises ValueError: If the property name is not a valid attribute.
- """
- if not hasattr(self, prop_name):
- raise ValueError("Unknown member variable given", prop_name)
- incr = test_equidistant(temp_array, 1.0e-3)
- if force_lookup or incr is False or len(temp_array) < temp_array_limit:
- # print("Performing perform_interpolation -> interpolate_lookup for", prop_name)
- result = interpolate_lookup(T, temp_array, prop_array)
- else:
- # print("Performing perform_interpolation -> interpolate_equidistant for", prop_name)
- result, sub_expr = interpolate_equidistant(
- T, float(temp_array[0]), float(incr), prop_array, prop_name)
- self.assignments.update({prop_name: sub_expr})
- setattr(self, prop_name, result)
+if __name__ == '__main__':
+ try:
+ # Valid case
+ Ti64 = Alloy([Ti, Al, V], [0.90, 0.06, 0.04], 1878, 1928)
+ print(f"Calculated Atomic Number: {0.90 * Ti.atomic_number + 0.06 * Al.atomic_number + 0.04 * V.atomic_number}")
+ print(f"Alloy Atomic Number: {Ti64.atomic_number}")
+ print(f"Initial Heat Conductivity: {Ti64.heat_conductivity}")
+ Ti64.heat_conductivity = 34
+ print(f"Updated Heat Conductivity: {Ti64.heat_conductivity}")
+ print(f"Boiling Temperature (Before Change): {Ti64.temperature_boil}")
+ Ti64.temperature_boil = 1000
+ print(f"Boiling Temperature (After Change): {Ti64.temperature_boil}")
+
+ # Invalid Composition
+ try:
+ invalid_alloy = Alloy([Ti, Al], [0.5, 0.5], 1878, 1928)
+ except ValueError as e:
+ print(f"Invalid Composition Test Passed: {e}")
+ # Empty Composition
+ try:
+ empty_composition_alloy = Alloy([Ti], [], 1878, 1928)
+ except ValueError as e:
+ print(f"Empty Composition Test Passed: {e}")
+
+ # Single Element Alloy
+ single_element_alloy = Alloy([Ti], [1.0], 1878, 1928)
+ print(f"Single Element Alloy Atomic Number: {single_element_alloy.atomic_number}")
+
+ # Invalid Property Assignment
+ try:
+ Ti64.heat_conductivity = "invalid_value"
+ except TypeError as e:
+ print(f"Invalid Property Assignment Test Passed: {e}")
+
+ # Boundary Values for Temperatures
+ boundary_alloy = Alloy([Ti, Al, V], [0.33, 0.33, 0.34], -273.15, 10000)
+ print(f"Boundary Temperatures: Solidus={boundary_alloy.temperature_solidus}, Liquidus={boundary_alloy.temperature_liquidus}")
+
+ # Properties Initialization
+ default_alloy = Alloy([Ti], [1.0], 1878, 1928)
+ print(f"Default Density: {default_alloy.density}")
+
+ except Exception as e:
+ print(f"An unexpected error occurred: {e}")
-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.atomic_number)
- print(Ti64.heat_conductivity)
- Ti64.heat_conductivity = 34
- print(Ti64.heat_conductivity)
- print(Ti64.temperature_boil)
- Ti64.temperature_boil = 1000
- print(Ti64.temperature_boil)
diff --git a/src/materials/alloys/SS316L/SS316L.py b/src/materials/alloys/SS316L/SS316L.py
index 06d3d726dbf7d2a82158d3ee49e76be7b99fec09..df963b13528362e1c72739af4ad5d180e8a592d1 100644
--- a/src/materials/alloys/SS316L/SS316L.py
+++ b/src/materials/alloys/SS316L/SS316L.py
@@ -1,28 +1,44 @@
-
-from pathlib import Path
-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.data_handler.data_handler import read_data, celsius_to_kelvin, thousand_times
-from src.materials.typedefs import VariableTypes
-from scipy.interpolate import interp1d
+from src.materials.interpolators import interpolate_property
-def create_SS316L(T: VariableTypes) -> Alloy:
+def create_SS316L(T: Union[float, sp.Symbol]) -> Alloy:
"""
Creates an Alloy instance for SS316L stainless steel with specific properties.
+ Args:
+ T (Union[float, sp.Symbol]): Temperature as a symbolic variable or numeric value.
+
Returns:
- Alloy: Initialized SS316L alloy with physical properties.
+ Tuple[Alloy, Tuple, Tuple, Tuple, Tuple]: Initialized SS316L alloy with physical properties and temperature arrays.
+
+ Notes:
+ - **Data Units**: All input data should be in SI units. If data is not in SI units, it must be converted. For example:
+ - **Temperature**: Convert Celsius to Kelvin using `celsius_to_kelvin` function.
+ - **Density**: Should be in kg/m³.
+ - **Heat Capacity**: Should be in J/(kg·K).
+ - **Heat Conductivity**: Should be in W/(m·K).
+ - **Temperature Array Consistency**: Ensure that all temperature arrays used for interpolation (`density_temp_array`, `heat_capacity_temp_array`, `heat_conductivity_temp_array`) have the same length.
+ In this implementation, `density_temp_array` is used as the reference array for all properties to ensure consistency.
+ - **Input Data Files**: The data files (`density_temperature.txt`, `heat_capacity_temperature.txt`, `heat_conductivity_temperature.txt`) must be located in the same directory as this script.
+ They should contain data in the units specified above.
+
+ 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.
"""
# Define the alloy with specific elemental composition and phase transition temperatures
SS316L = Alloy(
elements=[Fe, Cr, Mn, Ni],
composition=[0.708, 0.192, 0.018, 0.082], # Composition: 70.8% Fe, 19.2% Cr, 1.8% Mn, 8.2% Ni
temperature_solidus=1395.68, # Solidus temperature in Kelvin
- temperature_liquidus=1455.26 # Liquidus temperature in Kelvin
+ temperature_liquidus=1455.26, # Liquidus temperature in Kelvin
+ thermal_expansion_coefficient=1.7e-5 # in 1/K
)
# Determine the base directory
@@ -52,62 +68,18 @@ def create_SS316L(T: VariableTypes) -> Alloy:
heat_conductivity_temp_array = celsius_to_kelvin(heat_conductivity_temp_array)
- # Set the thermal expansion coefficient for SS316L
- SS316L.thermal_expansion_coefficient = 1.7e-5 # Thermal expansion coefficient in 1/K
-
# Perform interpolation for each property
- SS316L.interpolate_property(T, density_temp_array, density_array, 'density')
- # SS316L.interpolate_property(T, density_temp_array, density_array, 'density', tolerance=0.01) # Adaptive Interpolation
- SS316L.interpolate_property(T, heat_capacity_temp_array, heat_capacity_array, 'heat_capacity')
- SS316L.interpolate_property(T, heat_conductivity_temp_array, heat_conductivity_array, 'heat_conductivity')
-
- # Uncomment
- '''
- # Read temperatures for thermal diffusivity calculation
- thermal_diffusivity_temp_file_path = str(base_dir / 'thermal_diffusivity_temperature.txt')
- thermal_diffusivity_temp_array = np.asarray(read_temp(thermal_diffusivity_temp_file_path), dtype=float)
- '''
-
- # Create interpolation functions using scipy
- '''
- density_interp = interp1d(density_temp_array, density_array, kind='linear', fill_value='extrapolate')
- heat_capacity_interp = interp1d(heat_capacity_temp_array, heat_capacity_array, kind='linear', fill_value='extrapolate')
- heat_conductivity_interp = interp1d(heat_conductivity_temp_array, heat_conductivity_array, kind='linear', fill_value='extrapolate')
-
- # Evaluate numerical values
- heat_conductivity_values = heat_conductivity_interp(thermal_diffusivity_temp_array)
- density_values = density_interp(thermal_diffusivity_temp_array)
- heat_capacity_values = heat_capacity_interp(thermal_diffusivity_temp_array)
- '''
-
- # Create interpolation functions using numpy
- '''
- density_poly = np.polyfit(density_temp_array, density_array, deg=3) # Polynomial of degree 3
- density_func = np.poly1d(density_poly)
-
- heat_capacity_poly = np.polyfit(heat_capacity_temp_array, heat_capacity_array, deg=3) # Polynomial of degree 3
- heat_capacity_func = np.poly1d(heat_capacity_poly)
-
- heat_conductivity_poly = np.polyfit(heat_conductivity_temp_array, heat_conductivity_array, deg=3) # Polynomial of degree 3
- heat_conductivity_func = np.poly1d(heat_conductivity_poly)
- '''
-
- # Evaluate numerical values
- '''
- heat_conductivity_values = heat_conductivity_func(thermal_diffusivity_temp_array)
- density_values = density_func(thermal_diffusivity_temp_array)
- heat_capacity_values = heat_capacity_func(thermal_diffusivity_temp_array)
+ SS316L.density = interpolate_property(T, density_temp_array, density_array)
+ SS316L.heat_capacity = interpolate_property(T, heat_capacity_temp_array, heat_capacity_array)
+ SS316L.heat_conductivity = interpolate_property(T, heat_conductivity_temp_array, heat_conductivity_array)
# Calculate thermal diffusivity from thermal conductivity, density, and heat capacity
- thermal_diffusivity_values = thermal_diffusivity_by_heat_conductivity(heat_conductivity_values, density_values, heat_capacity_values)
-
- # Interpolate thermal diffusivity
- SS316L.interpolate_property(T, thermal_diffusivity_temp_array, np.asarray(thermal_diffusivity_values), "thermal_diffusivity")
- '''
-
- # Calculate thermal diffusivity from thermal conductivity, density, and heat capacity
- SS316L.thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(
- SS316L.heat_conductivity, SS316L.density, SS316L.heat_capacity)
+ 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)
+ ))
return SS316L
@@ -120,4 +92,28 @@ if __name__ == '__main__':
for i in range(len(alloy.composition)):
print(f"Element {alloy.elements[i]}: {alloy.composition[i]}")
- print(alloy.thermal_diffusivity)
+ print("\nTesting SS316L with symbolic temperature:")
+ for field in vars(alloy):
+ print(f"{field} = {alloy.__getattribute__(field)}")
+
+ # Test interpolate_property
+ print("\nTesting interpolate_property:")
+ test_temp = 1400.0 # Example temperature value
+
+ # Interpolate density, heat capacity, and heat conductivity
+ density_result = alloy.density.evalf(Temp, test_temp)
+ print(f"Interpolated density at T={test_temp} using evalf: {density_result}")
+
+ heat_capacity_result = alloy.heat_capacity.evalf(Temp, test_temp)
+ print(f"Interpolated heat capacity at T={test_temp} using evalf: {heat_capacity_result}")
+
+ heat_conductivity_result = alloy.heat_conductivity.evalf(Temp, test_temp)
+ print(f"Interpolated heat conductivity at T={test_temp} using evalf: {heat_conductivity_result}")
+
+ # 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 = thermal_diffusivity_by_heat_conductivity(
+ heat_conductivity, density, heat_capacity)
+ print(f"Calculated thermal diffusivity: {thermal_diffusivity}")
\ No newline at end of file
diff --git a/src/materials/alloys/Ti6Al4V.py b/src/materials/alloys/Ti6Al4V.py
index c49ed721b72cc8c87a6b4074b4edab27eda3239b..6ba95809333b53cd0da0f361f9de0438678f8218 100644
--- a/src/materials/alloys/Ti6Al4V.py
+++ b/src/materials/alloys/Ti6Al4V.py
@@ -1,68 +1,73 @@
-import numpy as np
import sympy as sp
+from typing import Union
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
+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.typedefs import VariableTypes
+from src.materials.typedefs import MaterialProperty
-def create_Ti6Al4V(T: VariableTypes) -> Alloy:
+def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
"""
Creates an Alloy instance for Ti6Al4V with specific properties.
Args:
- T (VariableTypes): Temperature as a symbolic variable or numeric value.
+ T (Union[float, sp.Symbol]): Temperature as a symbolic variable or numeric value.
Returns:
- Alloy: Initialized Ti6Al4V alloy with physical properties.
+ Alloy: Initialized Ti6Al4V alloy with physical properties including density, heat conductivity, and thermal diffusivity.
+
+ Notes:
+ - **Data Units**: All input data should be in SI units. If data is not in SI units, it must be converted. For example:
+ - **Temperature**: The temperature used in the function is expected to be in Kelvin. If you have temperature data in Celsius, convert it to Kelvin before using it.
+ - **Density**: Should be in kg/m³.
+ - **Heat Capacity**: Should be in J/(kg·K).
+ - **Heat Conductivity**: Should be in W/(m·K).
+ - **Temperature Array Consistency**: Ensure that all temperature arrays used for interpolation (`density_temp_array`, `heat_capacity_temp_array`, `heat_conductivity_temp_array`) have the same length to ensure consistency in property evaluation.
+ The current implementation uses the temperature range from the `solidification_interval` method for all interpolations to maintain consistency.
+ - **Solidus and Liquidus**: The solidus and liquidus temperatures define the temperature range over which the alloy is solid and liquid, respectively. These temperatures are used to determine the range for property interpolation.
+ - **Thermal Diffusivity Calculation**: The thermal diffusivity is computed using the formula:
+
+ Thermal Diffusivity = Heat Conductivity / Density * Heat Capacity
+
+ where:
+ - **Heat Conductivity**: Interpolated value based on temperature.
+ - **Density**: Computed based on thermal expansion.
+ - **Heat Capacity**: Provided directly.
+
+ 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.
+ For Ti6Al4V, the alloy properties are computed based on the provided temperature and constants, ensuring accurate representation of the material's behavior.
+
"""
- # Define alloy with specific elemental composition and phase transition temperatures
+ # Define the Ti6Al4V alloy with its elemental composition and phase transition temperatures
Ti6Al4V = Alloy(
elements=[Ti, Al, V],
- composition=[0.90, 0.06, 0.04], # 90% Ti, 6% Al, 4% V
- temperature_solidus=1878.0, # Temperature solidus in Kelvin
- temperature_liquidus=1928.0 # Temperature liquidus in Kelvin
+ composition=[0.90, 0.06, 0.04], # 90% Ti, 6% Al, 4% V
+ temperature_solidus=1878.0, # Temperature solidus in Kelvin
+ temperature_liquidus=1928.0, # Temperature liquidus in Kelvin
+ thermal_expansion_coefficient=8.6e-6, # 1/K [Source: MatWeb]
+ heat_capacity=526, # J/(kg·K) [Source: MatWeb]
+ latent_heat_of_fusion=290000.0, # J/kg [Source: MatWeb]
+ latent_heat_of_vaporization=8.86e6 # J/kg [Source: MatWeb]
)
- # Set properties specific to Ti6Al4V
- # Thermal expansion coefficient (1/K) [Source: MatWeb]
- Ti6Al4V.thermal_expansion_coefficient = 8.6e-6 # 1/K
-
- # Define density using models and interpolations for different temperature ranges
- # Density (kg/m³) [Source: MatWeb]
- # Ti6Al4V.density = 4430.0 # kg/m³
-
- # Using interpolate_equidistant for density
- T_array = np.linspace(Ti6Al4V.temperature_solidus, Ti6Al4V.temperature_liquidus, 7)
- # density_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
- density_array = [density_by_thermal_expansion(
- T_i, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient) for T_i in T_array]
- # Perform equidistant interpolation for the density at temperature T
- density_result, sub_expr_density = interpolate_equidistant(
- T, Ti6Al4V.temperature_solidus, 8, density_array, 'density')
- # Update the alloy with the interpolated density value
- Ti6Al4V.density = density_result
- # Store the sub-expressions for symbolic evaluation and further processing
- Ti6Al4V.assignments.update({"density": sub_expr_density})
-
- # Specific heat capacity (J/(kg·K)) [Source: MatWeb]
- Ti6Al4V.heat_capacity = 526 # J/(kg·K)
-
- # Thermal conductivity (W/(m·K)) [Source: MatWeb]
- Ti6Al4V.heat_conductivity = interpolate_lookup(
- T, Ti6Al4V.solidification_interval(), [6.7, 32.0]) # W/(m·K)
-
- # Latent heat of fusion (J/kg) [Source: MatWeb]
- Ti6Al4V.latent_heat_of_fusion = 290000.0 # J/kg
+ # Compute density based on thermal expansion and other constants
+ Ti6Al4V.density = MaterialProperty(density_by_thermal_expansion(
+ T, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient))
- # Latent heat of vaporization (J/kg) [Source: MatWeb]
- Ti6Al4V.latent_heat_of_vaporization = 8.86e6 # J/kg
+ # Interpolate heat conductivity based on temperature
+ Ti6Al4V.heat_conductivity = interpolate_property(
+ T, Ti6Al4V.solidification_interval(), [6.7, 32.0]) # W/(m·K)
- # Calculate thermal diffusivity from conductivity, density, and heat capacity
- Ti6Al4V.thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(
- Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
+ # Calculate thermal diffusivity using conductivity, density, and heat capacity
+ Ti6Al4V.thermal_diffusivity = interpolate_property(T, Ti6Al4V.solidification_interval(),
+ thermal_diffusivity_by_heat_conductivity(
+ Ti6Al4V.heat_conductivity.evalf(T, Ti6Al4V.solidification_interval()),
+ Ti6Al4V.density.expr,
+ Ti6Al4V.heat_capacity
+ ))
return Ti6Al4V
@@ -84,16 +89,13 @@ if __name__ == '__main__':
# Test interpolate_equidistant
print("\nTesting interpolate_equidistant:")
test_temp = 1850.0 # Example temperature value
- result, _ = 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 = interpolate_equidistant(
- Temp, 1820, 20.0, [700, 800, 900, 1000, 1100, 1200], 'heat_capacity')
- print(f"Assignments: {assignments}")
+ result = interpolate_equidistant(
+ Temp, 1820, 20.0, [700, 800, 900, 1000, 1100, 1200])
+ print(f"Interpolated heat capacity at T={test_temp} using interpolate_equidistant: {result.expr}")
+ print(f"Assignments: {result.assignments}")
- result_density, _ = interpolate_equidistant(
- test_temp, 1878, 8, [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], 'density')
+ result_density = interpolate_equidistant(
+ test_temp, 1878, 8, [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
print(f"Interpolated density at T={test_temp} using interpolate_equidistant: {result_density}")
# Test interpolate_lookup
@@ -111,8 +113,3 @@ if __name__ == '__main__':
thermal_diffusivity = thermal_diffusivity_by_heat_conductivity(
heat_conductivity, density, heat_capacity)
print(f"Calculated thermal diffusivity: {thermal_diffusivity}")
-
- # Uncomment below lines if plotting is needed
- # 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/assignment_converter.py b/src/materials/assignment_converter.py
index f8f0a304f5ee7db4ae5975e29335f210edb5af3b..043634c55363e556ef89c06feb88073a97b52bd1 100644
--- a/src/materials/assignment_converter.py
+++ b/src/materials/assignment_converter.py
@@ -1,10 +1,13 @@
import numpy as np
+import sympy as sp
import pystencils as ps
-from pystencils.types.quick import *
+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
-def type_mapping(type_str: str):
+def type_mapping(type_str: str) -> Union[np.dtype, Arr]:
"""
Maps a string representation of a type to a corresponding numpy or pystencils data type.
@@ -25,36 +28,71 @@ def type_mapping(type_str: str):
return np.dtype(type_str) # Default to numpy data type
-def assignment_converter(assignments_items):
+def assignment_converter(assignment_list: List[Assignment]) \
+ -> Tuple[List[ps.Assignment], Dict[sp.Symbol, ps.TypedSymbol]]:
"""
- Converts the assignment data from the Alloy class to pystencils-compatible format.
+ Converts a list of assignments from the Alloy class to pystencils-compatible format.
Parameters:
- assignments_items (Dict[str, List[Assignment]]): A dictionary where each key is a label (string) and each value is a list of `Assignment` objects.
- Each `Assignment` object should have the following attributes:
+ assignment_list (List[sp.Assignment]): A list of `Assignment` objects where each `Assignment` object should have:
- lhs (sp.Symbol): The left-hand side of the assignment, which is a symbolic variable. It should have a `name` and a `type` attribute.
- rhs (Union[tuple, sp.Expr]): The right-hand side of the assignment, which can be a tuple or a sympy expression.
- lhs_type (str): A string representing the type of the left-hand side, indicating the data type (e.g., "double[]", "float[]").
Returns:
- Tuple[List[ps.Assignment], Dict[Assignment, ps.TypedSymbol]]:
+ Tuple[List[ps.Assignment], Dict[sp.Symbol, ps.TypedSymbol]]:
- A list of `pystencils.Assignment` objects, each representing an assignment in the pystencils format.
- This list includes assignments where the `rhs` is directly assigned or cast to the appropriate type.
- - A dictionary mapping each original `Assignment` object to its corresponding `pystencils.TypedSymbol`.
- The dictionary maps `sp.Symbol` objects to `pystencils.TypedSymbol` instances, allowing for easy retrieval of types and assignments.
+ This list includes assignments where the `rhs` is directly assigned or cast to the appropriate type.
+ - A dictionary mapping each original `sp.Symbol` object (from the lhs) to its corresponding `pystencils.TypedSymbol`.
+ The dictionary allows for easy retrieval of types and assignments.
"""
subexp = []
subs = {}
- # Iterate over the dictionary items
- for label, assignment_list in assignments_items:
- print('assignment_list: ', assignment_list)
- for assignment in assignment_list:
- print('assignment in assignment_list: ', assignment)
- ps_type = type_mapping(assignment.lhs_type)
- data_symbol = ps.TypedSymbol(assignment.lhs.name, ps_type)
- if isinstance(ps_type, Arr):
- subexp.append(ps.Assignment(data_symbol, assignment.rhs))
- else: # needed for the sp.floor in interpolate_equidistant to cast it to an int from double
- subexp.append(ps.Assignment(data_symbol, CastFunc(assignment.rhs, ps_type)))
- subs[assignment.lhs] = data_symbol
+
+ for assignment in assignment_list:
+ ps_type = type_mapping(assignment.lhs_type)
+ data_symbol = ps.TypedSymbol(assignment.lhs.name, ps_type)
+
+ if isinstance(ps_type, Arr):
+ subexp.append(ps.Assignment(data_symbol, assignment.rhs))
+ else: # Cast rhs to the appropriate type when necessary
+ subexp.append(ps.Assignment(data_symbol, CastFunc(assignment.rhs, ps_type)))
+
+ subs[assignment.lhs] = data_symbol
+
return subexp, subs
+
+
+if __name__ == '__main__':
+ # Define test symbols and assignments
+ x = sp.Symbol('x')
+ y = sp.Symbol('y')
+ z = sp.Symbol('z')
+
+ # Example assignments
+ assignments = [
+ Assignment(lhs=x, rhs=sp.Symbol('value_x'), lhs_type='double[]'),
+ Assignment(lhs=y, rhs=sp.Symbol('value_y'), lhs_type='float[]'),
+ Assignment(lhs=z, rhs=sp.Symbol('value_z'), lhs_type='int')
+ ]
+
+ # 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
+
+ # Test assignment_converter function
+ print("\nTesting assignment_converter function:")
+ assignments_converted, symbol_map = assignment_converter(assignments)
+
+ for conv in assignments_converted:
+ print(conv)
+
+ print("\nSymbol to TypedSymbol mapping:")
+ for sym, typed_sym in symbol_map.items():
+ print(f"{sym}: {typed_sym}")
+
+ # Validate specific conditions for tests
+ assert isinstance(assignments_converted[0], ps.Assignment), "The converted assignment should be of type pystencils.Assignment"
+ assert isinstance(symbol_map[x], ps.TypedSymbol), "The symbol map should map symbols to pystencils.TypedSymbol"
diff --git a/src/materials/data_handler/data_handler.py b/src/materials/data_handler/data_handler.py
index b582b3c043100df0e65bb84c6adbf2aae6e364bc..1b4844f3b42e846a30c201a2e71ce47a4306d8e0 100644
--- a/src/materials/data_handler/data_handler.py
+++ b/src/materials/data_handler/data_handler.py
@@ -1,13 +1,11 @@
-import os
import numpy as np
from typing import Union, Tuple
-from scipy.interpolate import interp1d, CubicSpline, InterpolatedUnivariateSpline, Rbf
from src.materials.data_handler.create_test_files import create_test_files
def print_results(file_path: str, temperatures: np.ndarray, material_property: np.ndarray) -> None:
"""
- Prints the results.
+ Prints the results of the test.
Parameters:
file_path (str): The path to the data file.
@@ -20,107 +18,6 @@ def print_results(file_path: str, temperatures: np.ndarray, material_property: n
print("-" * 40)
-def interpolate_missing_values(temp: np.ndarray, prop: np.ndarray, target_temps: np.ndarray) \
- -> Tuple[np.ndarray, np.ndarray]:
- """
- Interpolates property values based on existing temperature data and target temperatures.
-
- Parameters:
- temp (np.ndarray): Array of existing temperature values.
- prop (np.ndarray): Array of property values corresponding to the temperatures.
- target_temps (np.ndarray): Array of target temperatures for interpolation.
-
- Returns:
- Tuple[np.ndarray, np.ndarray]: Target temperatures and interpolated property values.
- """
- # Perform linear interpolation using numpy.interp
- # interpolated_props = np.interp(target_temps, temp, prop, left=np.nan, right=np.nan)
- # return target_temps, interpolated_props
-
- # Perform linear interpolation using scipy.interpolate.interp1d
- interp_func = interp1d(temp, prop, kind='quadratic', fill_value='extrapolate')
- # interp_func = CubicSpline(temp, prop)
- # interp_func = InterpolatedUnivariateSpline(temp, prop)
- interpolated_props = interp_func(target_temps)
- return target_temps, interpolated_props
- # rbf = Rbf(temp, prop, function='multiquadric')
- # return target_temps, rbf(target_temps)
-
-
-def clean_temperature_data(temp: np.ndarray, prop: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
- """
- Cleans the temperature data by ensuring equally spaced values based on the most common temperature difference.
- Interpolates missing temperatures if necessary.
-
- Parameters:
- temp (np.ndarray): Array of temperature values.
- prop (np.ndarray): Array of property values corresponding to the temperatures.
-
- Returns:
- Tuple[np.ndarray, np.ndarray]: Cleaned temperatures and properties.
- """
- if len(temp) < 2:
- return temp, prop
-
- min_temp = np.min(temp)
- max_temp = np.max(temp)
-
- # Check and flip data if in descending order
- if np.all(np.diff(temp) < 0):
- temp = np.flip(temp)
- prop = np.flip(prop)
- print("Data was in descending order. Flipped to ascending.")
-
- diffs = np.diff(temp)
- print(f"Differences between consecutive temperatures: {diffs}")
-
- # Calculate the most common difference
- try:
- # common_diff = np.median(diffs) # Assuming the most common diff is the desired spacing
- # hist, bin_edges = np.histogram(diffs, bins='auto')
- # common_diff = bin_edges[np.argmax(hist)] # Most frequent temperature difference
- unique_diffs, counts = np.unique(diffs, return_counts=True)
- common_diff = unique_diffs[np.argmax(counts)]
- if np.max(counts) < 2:
- print("No repeated temperature differences found. Returning original data.")
- return temp, prop
- print(f"Most common temperature difference: {common_diff}")
- except Exception as e:
- print(f"Error calculating most common difference: {e}")
- return temp, prop
-
- if common_diff <= 0:
- raise ValueError("Temperature difference must be positive.")
-
- # start_temp = min_temp - (min_temp % common_diff) # Add boundary points by extrapolating out of bounds if required
- # expected_temps = np.arange(start_temp, max_temp + common_diff, common_diff) # Add boundary points by extrapolating out of bounds if required
-
- start_temp = min_temp + common_diff - (min_temp % common_diff) if min_temp % common_diff != 0 else min_temp # Cut boundary points if not equidistant
- end_temp = max_temp - (max_temp % common_diff) if max_temp % common_diff != 0 else max_temp # Cut boundary points if not equidistant
- expected_temps = np.arange(start_temp, end_temp + common_diff, common_diff)
-
- expected_temps = np.unique(expected_temps) # Remove duplicates
- print(f"Expected temperatures: {expected_temps}")
-
- # Handle case where expected_temps is empty
- if len(expected_temps) == 0:
- print("No expected temperatures generated, falling back to original temperatures.")
- return temp, prop
-
- # Interpolate properties for the expected temperatures
- cleaned_temp, cleaned_prop = interpolate_missing_values(temp, prop, expected_temps)
-
- print(f"Cleaned temperatures: {cleaned_temp}")
- print(f"Cleaned properties: {cleaned_prop}")
-
- # Ensure the final output is in ascending order
- sorted_indices = np.argsort(cleaned_temp)
- cleaned_temp = cleaned_temp[sorted_indices]
- cleaned_prop = cleaned_prop[sorted_indices]
-
- return cleaned_temp, cleaned_prop
-
-
def read_data(file_path: str, header: bool = True) -> Tuple[np.ndarray, np.ndarray]:
"""
Reads temperature and property data from a file and returns both the original and cleaned data.
@@ -161,16 +58,7 @@ def read_data(file_path: str, header: bool = True) -> Tuple[np.ndarray, np.ndarr
duplicate_rows = [str(idx + 1) for idx, value in enumerate(temp) if value in duplicates]
raise ValueError(f"Duplicate temperature entries found in rows: {', '.join(duplicate_rows)}")
- print(f"Original temperatures: {temp}")
- print(f"Original properties: {prop}")
-
- # Check if the data is already equidistant
- if np.all(np.diff(temp) == np.diff(temp)[0]):
- print("Data is already equidistant. No cleaning needed.")
- return temp, prop
-
- cleaned_temp, cleaned_prop = clean_temperature_data(temp, prop)
- return cleaned_temp, cleaned_prop
+ return temp, prop
except Exception as e:
print(f"Error reading data from {file_path}: {e}")
@@ -253,19 +141,21 @@ def run_tests() -> None:
"""
test_cases = [
{"name": "Normal Case", "file_path": "test_files/test_data_normal.txt"},
- # {"name": "Normal Case 1", "file_path": "test_files/test_data_normal1.txt"},
- # {"name": "Case with Variable Spacing", "file_path": "test_files/test_data_variable_spacing.txt"},
- # {"name": "Case with Sparse Data", "file_path": "test_files/test_data_sparse.txt"},
- # {"name": "Case with Missing Values", "file_path": "test_files/test_data_missing_values.txt"},
- # {"name": "Case with Descending Data", "file_path": "test_files/test_data_descending.txt"},
- # {"name": "Case with Duplicates", "file_path": "test_files/test_data_duplicates.txt"},
- # {"name": "Empty File", "file_path": "test_files/test_data_empty.txt"},
- # {"name": "File with Only Header", "file_path": "test_files/test_data_header.txt"},
- # {"name": "Invalid Data Entries", "file_path": "test_files/test_data_invalid_data.txt"},
- # {"name": "Mismatched Columns", "file_path": "test_files/test_data_mismatched_columns.txt"},
- # {"name": "Missing Data Entries", "file_path": "test_files/test_data_missing_data.txt"},
- # {"name": "Ascending density_temperature values", "file_path": "test_files/test_data_density_temperature_ascending.txt"},
- # {"name": "Descending density_temperature values", "file_path": "test_files/test_data_density_temperature_descending.txt"}
+ {"name": "Normal Case 1", "file_path": "test_files/test_data_normal1.txt"},
+ {"name": "Case with Variable Spacing", "file_path": "test_files/test_data_variable_spacing.txt"},
+ {"name": "Case with Sparse Data", "file_path": "test_files/test_data_sparse.txt"},
+ {"name": "Case with Missing Values", "file_path": "test_files/test_data_missing_values.txt"},
+ {"name": "Case with Descending Data", "file_path": "test_files/test_data_descending.txt"},
+ {"name": "Case with Duplicates", "file_path": "test_files/test_data_duplicates.txt"},
+ {"name": "Empty File", "file_path": "test_files/test_data_empty.txt"},
+ {"name": "File with Only Header", "file_path": "test_files/test_data_header.txt"},
+ {"name": "Invalid Data Entries", "file_path": "test_files/test_data_invalid_data.txt"},
+ {"name": "Mismatched Columns", "file_path": "test_files/test_data_mismatched_columns.txt"},
+ {"name": "Missing Data Entries", "file_path": "test_files/test_data_missing_data.txt"},
+ {"name": "Ascending density_temperature values",
+ "file_path": "test_files/test_data_density_temperature_ascending.txt"},
+ {"name": "Descending density_temperature values",
+ "file_path": "test_files/test_data_density_temperature_descending.txt"}
]
for case in test_cases:
@@ -296,5 +186,5 @@ def run_tests() -> None:
if __name__ == "__main__":
- create_test_files()
+ create_test_files() # Ensure test files are created before running tests
run_tests()
diff --git a/src/materials/interpolators.py b/src/materials/interpolators.py
index a308b8fb9ca8e34e77f03b75e895773f39cec04e..1a8b7341e26727a02e8a81eab79915f5b17e52f2 100644
--- a/src/materials/interpolators.py
+++ b/src/materials/interpolators.py
@@ -1,15 +1,16 @@
import numpy as np
import sympy as sp
-from typing import List, Tuple, Dict
-from src.materials.typedefs import *
-# from src.materials.data_handler.data_handler import test_equidistant
+from typing import Union
+from src.materials.typedefs import Assignment, ArrayTypes, MaterialProperty
+COUNT = 0
-def interpolate_equidistant(T: VariableTypes,
- T_base: float,
- T_incr: float,
- v_array: ValueTypes,
- label: str) -> Tuple[PropertyTypes, List[Assignment]]:
+
+def interpolate_equidistant(
+ T: Union[float, sp.Symbol],
+ T_base: float,
+ T_incr: float,
+ v_array: ArrayTypes) -> Union[float, MaterialProperty]:
"""
Perform equidistant interpolation for symbolic or numeric temperature values.
@@ -17,73 +18,75 @@ def interpolate_equidistant(T: VariableTypes,
:param T_base: Base temperature for the interpolation.
:param T_incr: Increment in temperature for each step in the interpolation array.
:param v_array: Array of values corresponding to the temperatures.
- :param label: Label for the symbolic representation.
- :return: Interpolated value and list of assignments (if symbolic).
+ :return: Interpolated value as a float or MaterialProperty object.
:raises ValueError: If the temperature type is unsupported.
"""
if T_incr < 0:
T_incr *= -1
- T_base = T_base - T_incr * (len(v_array) - 1)
+ T_base -= T_incr * (len(v_array) - 1)
v_array = np.flip(v_array)
if isinstance(T, sp.Symbol):
- # print("Performing interpolate_equidistant (symbolic)")
- sym_data = sp.IndexedBase(label, shape=len(v_array))
+ global COUNT
+ label = '_' + str(COUNT).zfill(3)
+ COUNT += 1
+
+ sym_data = sp.IndexedBase(label + "_var", shape=len(v_array))
sym_idx = sp.symbols(label + '_idx', cls=sp.Idx)
- sym_dec = sp.Symbol(label + "_pos_dec")
+ sym_dec = sp.Symbol(label + "_dec")
T_base = sp.Float(T_base)
T_incr = sp.Float(T_incr)
pos = (T - T_base) / T_incr
- pos_dec = pos - sp.floor(pos) # sp.frac(pos) not supported by pystencils yet!
+ 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),
- ((1 - sym_dec) * sym_data[sym_idx] + sym_dec * sym_data[sym_idx + 1], True))
+ ((1 - sym_dec) * sym_data[sym_idx] + sym_dec * sym_data[sym_idx + 1], True)
+ )
sub_expressions = [
Assignment(sym_data, tuple(sp.Float(float(v)) for v in v_array), "double[]"),
Assignment(sym_idx, pos, "int"),
Assignment(sym_dec, pos_dec, "double")
]
- return result, sub_expressions
+ return MaterialProperty(result, sub_expressions)
elif isinstance(T, float):
- # print("Performing interpolate_equidistant (numeric)")
v_array = np.asarray(v_array)
n = len(v_array)
min_temp = T_base
max_temp = T_base + (n - 1) * T_incr
if T < min_temp:
- # print(f"T ({T}) is below minimum temperature ({min_temp}), returning first value: {v_array[0]}.")
- return float(v_array[0]), []
+ return float(v_array[0])
elif T > max_temp:
- # print(f"T ({T}) is above maximum temperature ({max_temp}), returning last value: {v_array[-1]}.")
- return float(v_array[-1]), []
- pos = (T - T_base) / T_incr
+ return float(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]
- print(f"Interpolated value at T ({T}): {value}")
- return float(value), []
+ return float(value)
else:
raise ValueError(f"Unsupported type for T: {type(T)}")
-def interpolate_lookup(T: VariableTypes, T_array: ValueTypes, v_array: ValueTypes) -> PropertyTypes:
+def interpolate_lookup(
+ T: Union[float, sp.Symbol],
+ T_array: ArrayTypes,
+ v_array: ArrayTypes) -> Union[float, MaterialProperty]:
"""
Perform lookup-based interpolation for symbolic or numeric temperature values.
:param T: Temperature, either symbolic (sp.Symbol) or numeric (float).
:param T_array: Array of temperatures for lookup.
:param v_array: Array of values corresponding to the temperatures.
- :return: Interpolated value.
+ :return: Interpolated value as a float or MaterialProperty object.
:raises ValueError: If T_array and v_array lengths do not match or if T is of unsupported type.
"""
if len(T_array) != len(v_array):
@@ -94,9 +97,8 @@ def interpolate_lookup(T: VariableTypes, T_array: ValueTypes, v_array: ValueType
v_array = np.flip(v_array)
if isinstance(T, sp.Symbol):
- # print("Performing interpolate_lookup (symbolic)")
T_array = [sp.Float(float(t)) for t in T_array]
- v_array = [sp.Float(float(v)) for v in v_array]
+ v_array = [v if isinstance(v, sp.Expr) else sp.Float(float(v)) for v in v_array]
conditions = [
(v_array[0], T < T_array[0]),
@@ -106,10 +108,9 @@ def interpolate_lookup(T: VariableTypes, T_array: ValueTypes, v_array: ValueType
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)
+ return MaterialProperty(sp.Piecewise(*conditions))
elif isinstance(T, float):
- # print("Performing interpolate_lookup (numeric)")
T_array = np.asarray(T_array)
v_array = np.asarray(v_array)
if T < T_array[0]:
@@ -119,26 +120,75 @@ def interpolate_lookup(T: VariableTypes, T_array: ValueTypes, v_array: ValueType
else:
return float(np.interp(T, T_array, v_array))
- raise ValueError(f"Invalid input for T: Union[float, sp.Symbol] got {type(T)}")
+ else:
+ raise ValueError(f"Invalid input for T: {type(T)}")
+
+
+def test_equidistant(temp: np.ndarray) -> float:
+ """
+ Tests if the temperature values are equidistant.
+
+ :param temp: Array of temperature values.
+ :return: The common difference if equidistant, otherwise 0.
+ """
+ temperature_diffs = np.diff(temp)
+ if np.allclose(temperature_diffs, temperature_diffs[0], atol=1e-12):
+ return float(temperature_diffs[0])
+ return 0
+
+
+def interpolate_property(
+ T: Union[float, sp.Symbol],
+ temp_array: ArrayTypes,
+ prop_array: ArrayTypes,
+ temp_array_limit=6,
+ force_lookup=False) -> Union[float, MaterialProperty]:
+ """
+ Perform interpolation based on the type of temperature array (equidistant or not).
+
+ :param T: Temperature, either symbolic (sp.Symbol) or numeric (float).
+ :param temp_array: Array of temperatures.
+ :param prop_array: Array of property values corresponding to the temperatures.
+ :param temp_array_limit: Minimum length of the temperature array to force equidistant interpolation.
+ :param force_lookup: Boolean to force lookup interpolation regardless of temperature array type.
+ :return: Interpolated value.
+ :raises ValueError: If T_array and v_array lengths do not match.
+ """
+ if len(temp_array) != len(prop_array):
+ raise ValueError("T_array and v_array must have the same length")
+
+ if temp_array[0] > temp_array[-1]:
+ temp_array = np.flip(temp_array)
+ prop_array = np.flip(prop_array)
+
+ incr = test_equidistant(np.asarray(temp_array))
+
+ if force_lookup or incr == 0 or len(temp_array) < temp_array_limit:
+ return interpolate_lookup(T, temp_array, prop_array)
+ else:
+ return interpolate_equidistant(
+ T, float(temp_array[0]), float(incr), prop_array)
if __name__ == '__main__':
# Example usage and consistency tests for the interpolation functions
Temp = sp.Symbol('Temp')
- density_temp_array1 = np.array([600.0, 500.0, 400.0, 300.0, 200.0])
- density_v_array1 = np.array([10.0, 20.0, 30.0, 40.0, 50.0])
+ density_temp_array1 = np.flip(np.array([600.0, 500.0, 400.0, 300.0, 200.0]))
+ density_v_array1 = np.flip(np.array([10.0, 20.0, 30.0, 40.0, 50.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]
- Temp_base = 600.0
- Temp_incr = -100.0
+ Temp_base = 200.0
+ Temp_incr = 100.0
# 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)
+ print("Symbolic interpolation result with interpolate_lookup (arrays):",
+ symbolic_result.expr if isinstance(symbolic_result, MaterialProperty) else symbolic_result)
symbolic_result_list = interpolate_lookup(Temp, density_temp_array2, density_v_array2)
- print("Symbolic interpolation result with interpolate_lookup (lists):", symbolic_result_list)
+ print("Symbolic interpolation result with interpolate_lookup (lists):",
+ symbolic_result_list.expr if isinstance(symbolic_result_list, MaterialProperty) else symbolic_result_list)
numeric_result = interpolate_lookup(1894.7, density_temp_array2, density_v_array2)
print("Numeric interpolation result with interpolate_lookup (arrays):", numeric_result)
@@ -154,22 +204,22 @@ if __name__ == '__main__':
# Interpolation Equidistant Tests
print("\nInterpolate Equidistant Tests:")
- symbolic_result_eq, ass1 = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array1, 'density')
- print("Symbolic interpolation result with interpolate_equidistant (numpy arrays):", symbolic_result_eq)
- print("Assignments:", ass1)
+ 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, ass2 = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array2, 'density')
- print("Symbolic interpolation result with interpolate_equidistant (lists):", symbolic_result_eq_list)
- print("Assignments:", ass2)
+ 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, 'density')
+ 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)
- numeric_result_eq_list, _ = interpolate_equidistant(1900.2, Temp_base, Temp_incr, density_v_array2, 'density')
+ 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("invalid", Temp_base, Temp_incr, density_v_array2, 'density')
+ error_result_eq = interpolate_equidistant("what's up", Temp_base, Temp_incr, density_v_array2)
except ValueError as e:
print("Error:", e)
@@ -177,8 +227,9 @@ if __name__ == '__main__':
print("\nConsistency Test between interpolate_lookup and interpolate_equidistant:")
test_temps = [150.0, 200.0, 350.0, 450.0, 550.0, 650.0]
for T in test_temps:
- lookup_result = interpolate_lookup(T, 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(T, 600.0, -100.0, np.array([10.0, 20.0, 30.0, 40.0, 50.0]), 'density')
+ lookup_result = interpolate_lookup(T, 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(T, 200.0, 100.0, np.array([10.0, 20.0, 30.0, 40.0, 50.0]))
print(f"Temperature {T}:")
print(f"interpolate_lookup result: {lookup_result}")
print(f"interpolate_equidistant result: {equidistant_result}")
@@ -189,23 +240,8 @@ if __name__ == '__main__':
print(f"Skipping comparison for symbolic or non-numeric result at T = {T}")
symbolic_lookup_result = interpolate_lookup(Temp, density_temp_array1, density_v_array1)
- symbolic_equidistant_result, assignments = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array1, 'density')
+ symbolic_equidistant_result = interpolate_equidistant(Temp, Temp_base, Temp_incr, density_v_array1)
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}")
-
- # Perform Interpolation Tests
- # print("\nPerform Interpolation Tests:")
- # T_test = 350.0
- # temp_array = np.array([500.0, 400.0, 300.0, 200.0])
- # prop_array = np.array([40.0, 30.0, 20.0, 10.0])
- # assignments = {}
-
- # result_numeric = interpolate_pro(T_test, temp_array, prop_array, 'property', assignments)
- # print(f"Perform interpolation (numeric) result for T={T_test}: {result_numeric}")
-
- # T_symbolic = sp.Symbol('T')
- # result_symbolic = perform_interpolation(T_symbolic, temp_array, prop_array, 'property', assignments)
- # print(f"Perform interpolation (symbolic) result: {result_symbolic}")
- # print(f"Assignments: {assignments}")
+ print(f"interpolate_lookup (symbolic): {symbolic_lookup_result.expr}")
+ print(f"interpolate_equidistant (symbolic): {symbolic_equidistant_result.expr}")
+ print(f"interpolate_equidistant (symbolic): {symbolic_equidistant_result.assignments}")
diff --git a/src/materials/models.py b/src/materials/models.py
index f342c0404bd0f650c2e284bcaa91ae79e8e0c8e0..8a453662e1dbb8dab6d23ddce7825f13c531182f 100644
--- a/src/materials/models.py
+++ b/src/materials/models.py
@@ -1,13 +1,14 @@
+import numpy as np
import sympy as sp
-from src.materials.typedefs import PropertyTypes, VariableTypes, ValueTypes
+from typing import Union
def density_by_thermal_expansion(
- T: VariableTypes,
- base_temperature: float,
- base_density: float,
- thermal_expansion_coefficient: PropertyTypes) \
- -> PropertyTypes:
+ temperature: Union[float, np.ndarray, sp.Expr],
+ temperature_base: float,
+ density_base: float,
+ thermal_expansion_coefficient: Union[float, np.ndarray, sp.Expr]) \
+ -> Union[float, np.ndarray, sp.Expr]:
"""
Calculate density based on thermal expansion using the formula:
rho(T) = rho_0 / (1 + tec * (T - T_0))^3
@@ -17,20 +18,22 @@ def density_by_thermal_expansion(
- rho_0 is the density at the base temperature T_0.
- tec is the thermal expansion coefficient.
- :param T: Temperature value. Can be a numeric value (float) or a symbolic expression (sp.Symbol).
- :param base_temperature: Reference temperature at which the base density is given, as a float.
- :param base_density: Density at the reference temperature, as a float.
- :param thermal_expansion_coefficient: Thermal expansion coefficient. Can be a numeric value (float) or a symbolic expression (sp.Expr).
- :return: Density at temperature T. Returns a float if all inputs are numeric; otherwise, returns a sympy expression.
+ :param temperature: Temperature value. Can be a numeric value (float), numpy array (np.ndarray), or a symbolic expression (sp.Expr).
+ :param temperature_base: Reference temperature at which the base density is given, as a float.
+ :param density_base: Density at the reference temperature, as a float.
+ :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.
"""
- return base_density / (1 + thermal_expansion_coefficient * (T - base_temperature)) ** 3
+ 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)
def thermal_diffusivity_by_heat_conductivity(
- heat_conductivity: PropertyTypes,
- density: PropertyTypes,
- heat_capacity: PropertyTypes) \
- -> PropertyTypes:
+ heat_conductivity: Union[float, np.ndarray, sp.Expr],
+ density: Union[float, np.ndarray, sp.Expr],
+ heat_capacity: Union[float, np.ndarray, sp.Expr]) \
+ -> Union[float, np.ndarray, sp.Expr]:
"""
Calculate thermal diffusivity using the formula:
alpha(T) = k(T) / (rho(T) * c_p(T))
@@ -41,13 +44,60 @@ def thermal_diffusivity_by_heat_conductivity(
- rho(T) is the density at temperature T.
- c_p(T) is the specific heat capacity at temperature T.
- :param heat_conductivity: Thermal conductivity. Can be a numeric value (float) or a symbolic expression (sp.Expr).
- :param density: Density. Can be a numeric value (float) or a symbolic expression (sp.Expr).
- :param heat_capacity: Specific heat capacity. Can be a numeric value (float) or a symbolic expression (sp.Expr).
- :return: Thermal diffusivity. Returns a float if all inputs are numeric; otherwise, returns a sympy expression.
+ :param heat_conductivity: Thermal conductivity. Can be a numeric value (float), numpy array (np.ndarray), or a symbolic expression (sp.Expr).
+ :param density: Density. Can be a numeric value (float), numpy array (np.ndarray), or a symbolic expression (sp.Expr).
+ :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.
"""
+ 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)
- # Optional: Simplify the result if any input is an expression (sp.Expr)
- # if isinstance(heat_conductivity, sp.Expr) or isinstance(density, sp.Expr) or isinstance(heat_capacity, sp.Expr):
- # result = sp.simplify(result)
return result
+
+
+if __name__ == '__main__':
+ # Test 1: Single temperature value for density calculation
+ T_0 = 800
+ T_1 = 1000
+ tec = 1e-6
+ rho = 8000
+ print(f"Test 1 - Single temperature value: {T_1}, Density: {density_by_thermal_expansion(T_1, T_0, rho, tec)}")
+ # Expected output: Density value at T_1 = 1000
+
+ # Test 2: Temperature array for density calculation with constant TEC
+ T_a = np.linspace(1000, 2000, 3)
+ tec_a = np.ones(T_a.shape) * tec
+ print(f"Test 2a - Temperature array with constant TEC: {T_a}, Density: {density_by_thermal_expansion(T_a, T_0, rho, tec)}")
+ # Expected output: Array of densities for temperatures in T_a
+
+ # Test 3: Temperature array for density calculation with array TEC
+ print(f"Test 2b - Temperature array with array TEC: {T_a}, Density: {density_by_thermal_expansion(T_a, T_0, rho, tec_a)}")
+ # Expected output: Array of densities with temperature-dependent TEC
+
+ # Test 4: Thermal diffusivity calculation with scalar values
+ k = 30
+ c_p = 600
+ print(f"Test 3 - Thermal diffusivity with scalar values: {thermal_diffusivity_by_heat_conductivity(k, rho, c_p)}")
+ # Expected output: Scalar value of thermal diffusivity
+
+ # Test 5: Thermal diffusivity calculation with array values for heat conductivity
+ k_a = np.linspace(30, 40, 3)
+ print(f"Test 4a - Thermal diffusivity with array of heat conductivity: {thermal_diffusivity_by_heat_conductivity(k_a, rho, c_p)}")
+ # Expected output: Array of thermal diffusivities
+
+ # Test 6: Thermal diffusivity with density calculated by thermal expansion
+ calculated_densities = density_by_thermal_expansion(T_a, T_0, rho, tec)
+ print(f"Test 4b - Thermal diffusivity with calculated densities: {thermal_diffusivity_by_heat_conductivity(k_a, calculated_densities, c_p)}")
+ # Expected output: Array of thermal diffusivities considering temperature-dependent density
+
+ # Test 7: Symbolic computation for density with sympy Symbol temperature
+ T_symbolic = sp.Symbol('T')
+ print(f"Test 5 - Symbolic density computation: {density_by_thermal_expansion(T_symbolic, T_0, rho, tec)}")
+ # Expected output: Sympy expression for density as a function of temperature
+
+ # Test 8: Symbolic computation for thermal diffusivity
+ k_symbolic = sp.Symbol('k')
+ 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
diff --git a/src/materials/typedefs.py b/src/materials/typedefs.py
index 0fcd155cbcb4d3f8d62644af394c240f9370618d..7ba5ff70446c44c87c2cbf31ca62c22eae29ea84 100644
--- a/src/materials/typedefs.py
+++ b/src/materials/typedefs.py
@@ -1,30 +1,29 @@
"""
typedefs.py
-This module defines custom types and type aliases used throughout the materials module, enhancing
-code readability and consistency. It includes data structures for assignments and mappings
-between symbolic types and data types used in computations.
+This module defines custom types, type aliases, and utility classes used throughout the materials module,
+enhancing code readability and consistency. It includes data structures for assignments and the evaluation
+of material properties as functions of symbolic variables and numerical values.
Classes:
Assignment: A dataclass representing an assignment operation with a left-hand side (lhs),
right-hand side (rhs), and the type of the lhs.
+ MaterialProperty: A dataclass that represents a material property, which can be evaluated as a function
+ of a symbolic variable (e.g., temperature) and potentially includes symbolic assignments.
Type Aliases:
- ValueTypes: A type alias for numerical data types, which can be either numpy arrays or lists of floats.
- PropertyTypes: A type alias for properties that can be either constant floats or sympy expressions.
- VariableTypes: A type alias for variables that can be either floats or sympy symbols.
- AssignmentType: A type alias for a dictionary where keys are strings representing labels and values are lists of
- `Assignment` objects associated with those labels.
+ ArrayTypes: A type alias for numerical data types, which can be either numpy arrays, lists, or tuples.
+ PropertyTypes: A type alias for properties that can be either constant floats or MaterialProperty instances.
Functions:
- type_mapping(type_str: str) -> Union[np.dtype, Arr]:
- Maps a string representation of a type to a corresponding numpy or pystencils data type.
+ type_mapping(type_str: str) -> Union[np.dtype, Arr]: Maps a string representation of a type to a corresponding
+ numpy or pystencils data type.
"""
import numpy as np
import sympy as sp
-from dataclasses import dataclass
-from typing import Dict, List, Union
+from dataclasses import dataclass, field
+from typing import List, Union, Tuple, get_args
@dataclass
@@ -43,10 +42,156 @@ class Assignment:
# Type aliases for various data types used in the project
-ValueTypes = Union[np.ndarray, List[float]] # Numerical values can be represented as numpy arrays or lists of floats
-PropertyTypes = Union[float, sp.Expr] # Property values can be either constant floats or sympy expressions
-VariableTypes = Union[float, sp.Symbol] # Variables can be either numerical (floats) or symbolic (sympy symbols)
+ArrayTypes = Union[np.ndarray, List, Tuple] # Numerical values can be represented as numpy arrays, lists, or tuples
-# Type alias for a dictionary of assignments
-# Dictionary where keys are labels (strings) and values are lists of `Assignment` objects
-AssignmentType = Dict[str, List[Assignment]]
+
+@dataclass
+class MaterialProperty:
+ """
+ Represents a material property that can be evaluated as a function of a symbolic variable (e.g., temperature).
+ It can include symbolic assignments that are resolved during the evaluation.
+
+ Attributes:
+ expr (sp.Expr): The symbolic expression representing the material property.
+ assignments (List[Assignment]): A list of Assignment instances representing any symbolic assignments that
+ need to be resolved during evaluation.
+
+ Methods:
+ evalf(symbol: sp.Symbol, temperature: Union[float, ArrayTypes]) -> Union[float, np.ndarray]:
+ Evaluates the material property for a given temperature, resolving any symbolic assignments.
+ """
+ expr: sp.Expr
+ assignments: List[Assignment] = field(default_factory=list)
+
+ def evalf(self, symbol: sp.Symbol, temperature: Union[float, ArrayTypes]) -> Union[float, np.ndarray]:
+ """
+ Evaluates the material property at specific temperature values.
+
+ Args:
+ symbol (sp.Symbol): The symbolic variable (e.g., temperature) to substitute in the expression.
+ temperature (Union[float, ArrayTypes]): The temperature(s) at which to evaluate the property.
+
+ Returns:
+ Union[float, np.ndarray]: The evaluated property value(s) at the given temperature(s).
+ """
+ # If the expression has no symbolic variables, return it as a constant float
+ 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 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
+ elif isinstance(temperature, np.generic):
+ temperature = float(temperature)
+
+ # Prepare substitutions for symbolic assignments
+ substitutions = [(symbol, temperature)]
+ for assignment in self.assignments:
+ if isinstance(assignment.rhs, sp.Expr):
+ # Evaluate the right-hand side expression at the given temperature
+ value = sp.N(assignment.rhs.subs(symbol, temperature))
+ value = int(value) if assignment.lhs_type == "int" else value
+ substitutions.append((assignment.lhs, value))
+ else:
+ substitutions.append((assignment.lhs, assignment.rhs))
+
+ # Evaluate the material property with the substitutions
+ result = sp.N(self.expr.subs(substitutions))
+ return float(result)
+
+
+# Type alias for properties that can be either constant floats or MaterialProperty instances
+PropertyTypes = Union[float, MaterialProperty]
+
+if __name__ == '__main__':
+ # Example usage and tests for MaterialProperty and Assignment classes
+
+ T = sp.Symbol('T')
+ v = sp.IndexedBase('v')
+ i = sp.Symbol('i', type=sp.Integer)
+
+ # Example 1: Constant material property
+ mp0 = MaterialProperty(sp.Float(405.))
+ mp0.assignments.append(Assignment(sp.Symbol('A'), (100, 200), 'int'))
+ print(mp0)
+ print(mp0.evalf(T, 100.))
+
+ # Example 2: Linear temperature-dependent property
+ mp1 = MaterialProperty(T * 100.)
+ print(mp1)
+ print(mp1.evalf(T, 100.))
+
+ # Example 3: Indexed base with symbolic assignments
+ mp2 = MaterialProperty(v[i])
+ mp2.assignments.append(Assignment(v, (3, 6, 9), 'float'))
+ mp2.assignments.append(Assignment(i, T / 100, 'int'))
+ print(mp2)
+ print(mp2.evalf(T, 97)) # Should evaluate with i=0 (since 97/100 is 0 when converted to int)
+ print(mp2.evalf(T, 107)) # Should evaluate with i=1 (since 107/100 is 1 when converted to int)
+
+ # Example 4: Evaluate an expression with an array of temperatures
+ rho = 1e-6 * T ** 3 * 4000
+ print(rho * np.array([10, 20, 30]))
+
+if __name__ == '__main__':
+ # Example usage and tests for MaterialProperty and Assignment classes
+
+ T = sp.Symbol('T')
+ v = sp.IndexedBase('v')
+ i = sp.Symbol('i', type=sp.Integer)
+
+ # Test 1: Constant material property
+ mp0 = MaterialProperty(sp.Float(405.))
+ mp0.assignments.append(Assignment(sp.Symbol('A'), (100, 200), 'int'))
+ print(f"Test 1 - Constant property, no symbolic variables: {mp0.evalf(T, 100.)}") # Expected output: 405.0
+
+ # Test 2: Linear temperature-dependent property
+ mp1 = MaterialProperty(T * 100.)
+ print(f"Test 2 - Linear temperature-dependent property: {mp1.evalf(T, 100.)}") # Expected output: 10000.0
+
+ # Test 3: Indexed base with symbolic assignments
+ mp2 = MaterialProperty(v[i])
+ mp2.assignments.append(Assignment(v, (3, 6, 9), 'float'))
+ mp2.assignments.append(Assignment(i, T / 100, 'int'))
+ print(f"Test 3a - Indexed base with i=0: {mp2.evalf(T, 97)}") # Expected output: 3 (i=0)
+ print(f"Test 3b - Indexed base with i=1: {mp2.evalf(T, 107)}") # Expected output: 6 (i=1)
+
+ # Test 4: Evaluate an expression with an array of temperatures
+ rho = 1e-6 * T ** 3 * 4000
+ temps = np.array([10, 20, 30])
+ print(f"Test 4 - Expression evaluated with array: {rho * temps}") # Expected output: array of evaluated values
+
+ # Additional Tests:
+
+ # Test 5: Non-linear temperature-dependent property
+ mp3 = MaterialProperty(T ** 2 + T * 50 + 25)
+ print(f"Test 5 - Non-linear temperature-dependent property: {mp3.evalf(T, 10)}") # Expected output: 625.0
+
+ # Test 6: Temperature array evaluation for non-linear expression
+ print(f"Test 6 - Non-linear property with temperature array: {mp3.evalf(T, temps)}")
+ # Expected output: array of evaluated values for T in [10, 20, 30]
+
+ # Test 7: Property with multiple symbolic assignments
+ mp4 = MaterialProperty(v[i] * T)
+ mp4.assignments.append(Assignment(v, (3, 6, 9), 'float'))
+ mp4.assignments.append(Assignment(i, T / 100, 'int'))
+ print(f"Test 7a - Property with multiple symbolic assignments at T=95: {mp4.evalf(T, 95)}") # Expected: 285.0
+ print(f"Test 7b - Property with multiple symbolic assignments at T=205: {mp4.evalf(T, 205)}") # Expected: 1230.0
+
+ # Test 8: Constant property evaluated with temperature array
+ mp5 = MaterialProperty(sp.Float(500.))
+ print(f"Test 8 - Constant property with temperature array: {mp5.evalf(T, temps)}")
+ # Expected output: array of 500s for each temperature
+
+ # Test 9: Handling numpy scalar input
+ scalar_temp = np.float64(150.0)
+ mp6 = MaterialProperty(T + 50)
+ print(f"Test 9 - Handling numpy scalar input: {mp6.evalf(T, scalar_temp)}") # Expected output: 200.0
+
+ # 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