working materials package using property.py

apps/tutorials/codegen/HeatEquationKernel.py

 import sys
 sys.path.append('/local/ca00xebo/repos/walberla')
-import numpy as np
 import sympy as sp
 import pystencils as ps
 from pystencils_walberla import CodeGeneration, generate_sweep
-from dataclasses import dataclass
-from pystencils.typing import PointerType
-import src.materials.alloys.Ti6Al4V as ti
-from src.materials.get_parameters import get_parameters
-from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
-
-'''
-def interpolate(x, xmin, xmax, ymin, ymax):
-    """
-    Linear interpolation between two points: y = y_1 + ((y2 - y1) / (x2 - x1)) * (x - x1)
-    """
-    return (ymin * (xmax - x) + ymax * (x - xmin)) / (xmax - xmin)
-
-
-def interpolate_property(T, base_temperature, base_values):  # TODO
-    """
-    Interpolates properties based on temperature for both symbolic and numeric cases.
-
-    :param T: Temperature value, either symbolic or numeric.
-    :param base_temperature: List or array of base temperatures.
-    :param base_values: List or array of property values corresponding to base temperatures.
-    """
-    # Ensure base_temperature and base_values are lists or arrays
-    if not isinstance(base_temperature, (list, np.ndarray)) or not isinstance(base_values, (list, np.ndarray)):
-        raise TypeError("base_temperature and base_values must be of type list or np.ndarray")
-    # Convert to lists if they are numpy arrays
-    if isinstance(base_temperature, np.ndarray):
-        base_temperature = base_temperature.tolist()
-    if isinstance(base_values, np.ndarray):
-        base_values = base_values.tolist()
-    # Check if base_temperature and base_values have the same length
-    if len(base_temperature) != len(base_values):
-        raise ValueError("base_temperature and base_values must have the same length")
-    # assert len(base_temperature) == len(base_values), "base_temperature and base_values must have the same length"
-
-    if isinstance(T, sp.Expr):  # Symbolic expression
-        base_temperature = [sp.Float(t) for t in base_temperature]
-        base_values = [sp.Float(val) for val in base_values]
-        conditions = [
-            (base_values[0], T < base_temperature[0]),
-            (base_values[-1], T >= base_temperature[-1])]
-        for i in range(len(base_temperature) - 1):
-            interp_expr = (base_values[i] +
-                           (base_values[i+1] - base_values[i]) / (base_temperature[i+1] - base_temperature[i]) *
-                           (T - base_temperature[i]))
-            conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
-        return sp.Piecewise(*conditions)
-    else:  # Numeric expression
-        return np.interp(T, base_temperature, base_values)
-
-
-def density_by_thermal_expansion(T, base_temperature, base_density, thermal_expansion_coefficient):
-    """
-    Calculate density based on thermal expansion: rho_0 / (1 + tec * (T - T_0))
-
-    :param T: Temperature value.
-    :param base_temperature: Reference temperature.
-    :param base_density: Density at reference temperature.
-    :param thermal_expansion_coefficient: Thermal expansion coefficient.
-    """
-    return base_density / (1 + thermal_expansion_coefficient * (T - base_temperature)) ** 3
-
-
-def thermal_diffusivity_func(heat_conductivity, density, specific_heat_capacity):
-    """
-    Calculate thermal diffusivity: k(T) / (rho(T) * c_p(T))
-
-    :param heat_conductivity: Thermal conductivity.
-    :param density: Density.
-    :param specific_heat_capacity: Specific heat capacity.
-    """
-    return heat_conductivity / (density * specific_heat_capacity)
-
-
-def density_lookup(T, temperature_array, density_array):
-    """
-    Look up density based on temperature using interpolation.
-
-    :param T: Temperature value.
-    :param temperature_array: Array of temperatures.
-    :param density_array: Array of densities.
-    """
-    if isinstance(T, sp.Expr):  # Symbolic expression
-        conditions = [
-            (density_array[0], T < temperature_array[0]),
-            (density_array[-1], T >= temperature_array[-1])
-        ]
-        for i in range(len(temperature_array) - 1):
-            interp_expr = (density_array[i] +
-                           (density_array[i+1] - density_array[i]) / (temperature_array[i+1] - temperature_array[i]) *
-                           (T - temperature_array[i]))
-            conditions.append((interp_expr, sp.And(T >= temperature_array[i], T < temperature_array[i + 1])))
-        return sp.Piecewise(*conditions)
-    else:  # Numeric expression
-        return np.interp(T, temperature_array, density_array)
-
-
-def get_arguments(T, temperatures, args):
-    """
-    Recursively get arguments for parameter functions.
-
-    :param T: Temperature value.
-    :param temperatures: Solidus and liquidus temperatures.
-    :param args: Arguments to be processed.
-    """
-    converted_args = []
-    for arg in args:
-        if isinstance(arg, (Solidus_Liquidus_Parameter, Functional_Parameter)):
-            converted_args.append(get_parameters(T, temperatures, arg))
-        else:
-            converted_args.append(arg)
-    return converted_args
-
-
-def get_parameters(T, temperatures, parameter):
-    """
-    Resolve parameters based on temperature and parameter type.
-
-    :param T: Temperature value.
-    :param temperatures: Solidus and liquidus temperatures.
-    :param parameter: Parameter to be resolved.
-    """
-    if isinstance(parameter, Solidus_Liquidus_Parameter):
-        return sp.Piecewise(
-            (parameter.sol, T <= temperatures.sol),
-            (parameter.liq, T >= temperatures.liq),
-            (interpolate(T, temperatures.sol, temperatures.liq, parameter.sol, parameter.liq), True)
-        )
-    elif isinstance(parameter, Functional_Parameter):
-        resolved_args = [get_parameters(T, temperatures, arg) for arg in parameter.args]
-        if parameter.func == 'thermal_expansion':
-            return density_by_thermal_expansion(T, *resolved_args)  # *get_parameters(T, temperatures, parameter.args)
-        elif parameter.func == 'interpolation':
-            return interpolate_property(T, *resolved_args)
-        elif parameter.func == 'interpolation_in_array':  # TODO
-            T_base, T_incr, data_array, data_symbol = resolved_args
-            idx = sp.floor((T - T_base) / T_incr)
-            mod = (T - T_base) - idx * T_incr  # mod replaces (T - T_base) % T_incr
-            return sp.Piecewise(
-                (data_array[0], T < T_base),
-                (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
-                (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True)
-            )
-        elif parameter.func == 'lookup':
-            return density_lookup(T, *resolved_args)
-        elif parameter.func == 'thermal_diffusivity':
-            return thermal_diffusivity_func(*resolved_args)
-        else:
-            raise ValueError("No known format in the dict: parameter = ", parameter)
-    else:
-        return parameter  # assume everything which is not a dict is a float, int, etc.
-
-
-# Material Data Class
-@dataclass
-class Ti6Al4V:                                          # 90% Ti, 6% Al, 4% V
-    latent_heat =                           290000      # m²/s²
-    temperature_solidus_1 =                 1876.0      # k         Ref [1]
-    temperature_solidus_2 =                 1880.0      # k         Ref [2]
-    temperature_liquidus =                  1928.0      # k         Ref [2]
-    density_room =                          4430.0      # kg/m³     Ref [5]
-    thermal_expansion_coefficient =         10e-6       # /k        Ref[6]
-    heat_capacity_solidus =                 700.0       # J/kg-k
-    heat_capacity_liquidus =                1126.0      # J/kg-k
-    heat_conductivity_solidus =             21.58       # W/m-k
-    heat_conductivity_liquidus =            33.60       # W/m-k
-    density_temp_array = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
-    density_data_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
-    density_temp_base =                     1878.0
-    density_temp_incr =                     8.0
-    # thermal_diffusivity_solidus = 727.7
-    # thermal_diffusivity_liquidus = 708.6
-    # electrical_resistivity = [1, 2, 5, 10, 20, 50, 100]
-    # electrical_resistivity_temperatures = [800, 900, 1200, 1500, 1700, 1900, 2000]
-    # Define temperature and density arrays
-    density_temp_array_lookup = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
-    density_data_array_lookup = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
-
-
-@dataclass
-class Constants:
-    temperature_room =                      298.15          # k
-    N_a =                                   6.022141e23     # /mol
-    u =                                     1.660538e-27    # kg
-
-
-# Instantiate the Ti6Al4V class
-ti6al4v = Ti6Al4V()
-
-Ti6Al4V_latent_heat = ti6al4v.latent_heat
-Ti6Al4V_thermal_expansion_coefficient = ti6al4v.thermal_expansion_coefficient
-
-
-@dataclass
-class Solidus_Liquidus_Parameter:
-    sol: float
-    liq: float
-
-
-@dataclass
-class Functional_Parameter:
-    func: str
-    args: list
-
-
-# Parameters for Ti6Al4V
-Ti6Al4V_temperatures = Solidus_Liquidus_Parameter(
-    sol=0.5 * (ti6al4v.temperature_solidus_1 + ti6al4v.temperature_solidus_2),
-    liq=ti6al4v.temperature_liquidus)
-
-Ti6Al4V_heat_conductivity = Solidus_Liquidus_Parameter(
-    sol=ti6al4v.heat_conductivity_solidus,
-    liq=ti6al4v.heat_conductivity_liquidus)
-
-Ti6Al4V_specific_heat_capacity = Solidus_Liquidus_Parameter(
-    sol=ti6al4v.heat_capacity_solidus,
-    liq=ti6al4v.heat_capacity_liquidus)
-
-Ti6Al4V_density_tec = Functional_Parameter(
-    func='thermal_expansion',
-    args=[Constants.temperature_room, ti6al4v.density_room, ti6al4v.thermal_expansion_coefficient])
-
-Ti6Al4V_density_int = Functional_Parameter(
-    func='interpolation',
-    args=[ti6al4v.density_temp_array, ti6al4v.density_data_array])
-
-Ti6Al4V_density_in_array = Functional_Parameter(
-    func='interpolation_in_array',  # TODO
-    args=[ti6al4v.density_temp_base, ti6al4v.density_temp_incr, ti6al4v.density_data_array,
-          sp.IndexedBase(ps.TypedSymbol('density_data_symbol', PointerType(sp.Basic('float32'))),  # PointerType(create_type(np.float32)) PointerType(base_type=float32) PointerType(sp.Basic('float32'))
-                         shape=len(ti6al4v.density_data_array))])
-
-Ti6Al4V_density_lookup = Functional_Parameter(
-    func='lookup',
-    args=[ti6al4v.density_temp_array_lookup, ti6al4v.density_data_array_lookup])
-
-Ti6Al4V_thermal_diffusivity = Functional_Parameter(
-    func='thermal_diffusivity',
-    args=[Ti6Al4V_heat_conductivity, Ti6Al4V_density_lookup, Ti6Al4V_specific_heat_capacity]
-    # solidus=ti6al4v.thermal_diffusivity_solidus,
-    # liquidus=ti6al4v.thermal_diffusivity_liquidus
-)
-'''
+from src.materials.alloys import Ti6Al4V
+from src.materials.models import thermal_diffusivity_by_heat_conductivity
+from pystencils.types.quick import *
 with CodeGeneration() as ctx:
     data_type = "float64" if ctx.double_accuracy else "float32"
 
@@ -260,18 +19,25 @@ with CodeGeneration() as ctx:
     discretize = ps.fd.Discretization2ndOrder(dx=dx, dt=dt)
     heat_pde_discretized = discretize(heat_pde)
     heat_pde_discretized = heat_pde_discretized.args[1] + heat_pde_discretized.args[0].simplify()
-
-    thermal_diffusivity_expr = get_parameters(u.center(), Ti6Al4V_parameters.Ti6Al4V_temperatures, Ti6Al4V_parameters.Ti6Al4V_thermal_diffusivity)
-
+    Ti6Al4V = Ti6Al4V.create_Ti6Al4V(u.center())
+    # density, density_subexpr = Ti6Al4V.density  # Alternative
+    subexp = []
+    for a in Ti6Al4V.assignments:
+        data_symbol = ps.TypedSymbol(a[0].label.name, Arr(Fp(64)))  # Alternative
+        subexp.append(ps.Assignment(data_symbol, a[1]))
+    # thermal_diffusivity_expr = thermal_diffusivity_by_heat_conductivity(Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
+    thermal_diffusivity_expr = thermal_diffusivity_by_heat_conductivity(Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)  # Alternative
+    subexp.append(ps.Assignment(thermal_diffusivity, thermal_diffusivity_expr))
     ac = ps.AssignmentCollection(
-        subexpressions=[
-            ps.Assignment(thermal_diffusivity, thermal_diffusivity_expr)
-        ],
+        subexpressions=subexp,
         main_assignments=[
             ps.Assignment(u_tmp.center(), heat_pde_discretized),
             ps.Assignment(thermal_diffusivity_out.center(), thermal_diffusivity)
         ])
-    ac = ps.simp.simplifications.add_subexpressions_for_divisions(ac)
+
+    # ac = ac.new_with_substitutions(subexp)
+
+    # 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)
 

src/materials/alloys/Alloy.py

+import numpy as np
+from sympy import Expr
+from typing import List, Union
+from dataclasses import dataclass, field, asdict
+from src.materials.elements import Element
+from src.materials.elements import Ti, Al, V
+
+
+# Alloy Base Class
+@dataclass
+class Alloy:
+    elements: List[Element]
+    composition: List[float]
+
+    atomic_number: float = field(init=False)
+    atomic_mass: float = field(init=False)
+    temperature_boil: float = field(init=False)
+
+    temperature_solidus: float
+    temperature_liquidus: float
+
+    latent_heat_of_fusion: Union[float, Expr] = None
+    latent_heat_of_vaporization: Union[float, Expr] = None
+    density: Union[float, Expr] = None
+    thermal_expansion_coefficient: Union[float, Expr] = None
+    heat_capacity: Union[float, Expr] = None
+    heat_conductivity: Union[float, Expr] = None
+    thermal_diffusivity: Union[float, Expr] = None
+
+    assignments: list[tuple, ...] = field(default_factory=list, init=False)
+
+    def interpolate(self, quantity: str) -> float:
+        comp = np.asarray(self.composition)
+        qnt = np.asarray([element.__getattribute__(quantity) for element in self.elements])
+        return float(np.sum(comp * qnt))
+
+    def __post_init__(self):
+        self.atomic_number = self.interpolate("atomic_number")
+        self.atomic_mass = self.interpolate("atomic_mass")
+        self.temperature_boil = self.interpolate("temperature_boil")
+
+    def solidification_interval(self) -> [float, float]:
+        return [self.temperature_solidus, self.temperature_liquidus]
+
+
+if __name__ == '__main__':
+    Ti64 = Alloy([Ti, Al, V], [0.90, 0.06, 0.04], 1878, 1928)
+    print(0.9*22+0.06*13+0.04*23, Ti64.interpolate("atomic_number"), Ti64.atomic_number)
+    print(Ti64.heat_conductivity)
+    Ti64.heat_conductivity = 34
+    print(Ti64.heat_conductivity)
+    print(Ti64.temperature_boil)
+    Ti64.temperature_boil = 1000
+    print(Ti64.temperature_boil)

src/materials/alloys/Ti6Al4V.py

+import sympy as sp
+from typing import Union
+from src.materials.alloys.Alloy import Alloy
+from src.materials.elements import Ti, Al, V
+from src.materials.constants import Constants
+from src.materials import interpolators
+from src.materials import models
+from src.materials.property import Property
+
+
+def create_Ti6Al4V(T: Union[float, sp.Symbol]) -> Alloy:
+    # 90% Ti, 6% Al, 4% V
+    Ti6Al4V = Alloy(
+        [Ti, Al, V],
+        [0.90, 0.06, 0.04],
+        1878.0,
+        1928.0
+    )
+
+    Ti6Al4V.latent_heat_of_fusion = 290000.0  # m²/s²
+    Ti6Al4V.density = 4430.0  # kg/m³     Ref [5]
+    Ti6Al4V.thermal_expansion_coefficient = 10e-6  # /K        Ref[6]
+
+    '''Ti6Al4V.heat_capacity = interpolators.interpolate_lookup(
+        T, Ti6Al4V.solidification_interval(), [700.0, 1126.0])  # J/kg-k'''
+
+    heat_capacity = interpolators.interpolate_equidistant(
+        T, 1820, 20, [700, 800, 900, 1000, 1100, 1200], 'heat_capacity')  # J/kg-k
+
+    Ti6Al4V.heat_capacity = heat_capacity.expression
+    Ti6Al4V.assignments += heat_capacity.assignments
+
+    Ti6Al4V.heat_conductivity = interpolators.interpolate_lookup(
+        T, Ti6Al4V.solidification_interval(), [21.58, 33.60])  # W/m-k
+
+    '''Ti6Al4V.density = models.density_by_thermal_expansion(
+        T, Constants.temperature_room, 4430, Ti6Al4V.thermal_expansion_coefficient)
+
+    Ti6Al4V.density = interpolators.interpolate_lookup(
+        T, [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0],
+        [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])'''
+
+    density = interpolators.interpolate_equidistant(
+        T, 1878.0, 5,
+        [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0], 'density')
+
+    Ti6Al4V.density = density.expression
+    Ti6Al4V.assignments += density.assignments
+
+    Ti6Al4V.thermal_diffusivity = models.thermal_diffusivity_by_heat_conductivity(
+        Ti6Al4V.heat_conductivity, Ti6Al4V.density, Ti6Al4V.heat_capacity)
+
+    return Ti6Al4V
+
+
+if __name__ == '__main__':
+    Temp = sp.Symbol('Temp')
+    alloy = create_Ti6Al4V(Temp)
+    print('Ti6Al4V')
+    for field in vars(alloy):
+        print(field, "=", alloy.__getattribute__(field))
+
+    print(models.thermal_diffusivity_by_heat_conductivity(500, 5000, 600))
+
+    print(alloy.density.subs(Temp, 1890))

src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py

-import numpy as np
-import sympy as sp
-import pystencils as ps
-from dataclasses import dataclass
-from pystencils.typing import PointerType
-from src.materials.constants import Constants
-from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
-
-
-# Material Data Class
-@dataclass
-class Ti6Al4V:                                          # 90% Ti, 6% Al, 4% V
-    latent_heat =                           290000      # m²/s²
-    temperature_solidus_1 =                 1876.0      # k         Ref [1]
-    temperature_solidus_2 =                 1880.0      # k         Ref [2]
-    temperature_liquidus =                  1928.0      # k         Ref [2]
-    density_room =                          4430.0      # kg/m³     Ref [5]
-    thermal_expansion_coefficient =         10e-6       # /k        Ref[6]
-    heat_capacity_solidus =                 700.0       # J/kg-k
-    heat_capacity_liquidus =                1126.0      # J/kg-k
-    heat_conductivity_solidus =             21.58       # W/m-k
-    heat_conductivity_liquidus =            33.60       # W/m-k
-    density_temp_array = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
-    density_data_array = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
-    density_temp_base =                     1878.0
-    density_temp_incr =                     8.0
-    thermal_diffusivity_solidus =           727.7
-    thermal_diffusivity_liquidus =          708.6
-    # electrical_resistivity = [1, 2, 5, 10, 20, 50, 100]
-    # electrical_resistivity_temperatures = [800, 900, 1200, 1500, 1700, 1900, 2000]
-    # Define temperature and density arrays
-    density_temp_array_lookup = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
-    density_data_array_lookup = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
-
-
-'''
-@dataclass
-class Constants:
-    temperature_room =                      298.15          # k
-    N_a =                                   6.022141e23     # /mol
-    u =                                     1.660538e-27    # kg
-'''
-
-# Instantiate the Ti6Al4V class
-ti6al4v = Ti6Al4V()
-
-Ti6Al4V_latent_heat = ti6al4v.latent_heat
-Ti6Al4V_thermal_expansion_coefficient = ti6al4v.thermal_expansion_coefficient
-
-'''
-@dataclass
-class Solidus_Liquidus_Parameter:
-    sol: float
-    liq: float
-
-
-@dataclass
-class Functional_Parameter:
-    func: str
-    args: list
-'''
-
-# Parameters for Ti6Al4V
-Ti6Al4V_temperatures = Solidus_Liquidus_Parameter(
-    sol=0.5 * (ti6al4v.temperature_solidus_1 + ti6al4v.temperature_solidus_2),
-    liq=ti6al4v.temperature_liquidus)
-
-Ti6Al4V_heat_conductivity = Solidus_Liquidus_Parameter(
-    sol=ti6al4v.heat_conductivity_solidus,
-    liq=ti6al4v.heat_conductivity_liquidus)
-
-Ti6Al4V_specific_heat_capacity = Solidus_Liquidus_Parameter(
-    sol=ti6al4v.heat_capacity_solidus,
-    liq=ti6al4v.heat_capacity_liquidus)
-
-Ti6Al4V_density_tec = Functional_Parameter(
-    func='thermal_expansion',
-    args=[Constants.temperature_room, ti6al4v.density_room, ti6al4v.thermal_expansion_coefficient])
-
-Ti6Al4V_density_int = Functional_Parameter(
-    func='interpolation',
-    args=[ti6al4v.density_temp_array, ti6al4v.density_data_array])
-
-Ti6Al4V_density_in_array = Functional_Parameter(
-    func='interpolation_in_array',  # TODO
-    args=[ti6al4v.density_temp_base, ti6al4v.density_temp_incr, ti6al4v.density_data_array,
-          sp.IndexedBase(ps.TypedSymbol('density_data_symbol', PointerType(sp.Basic('float32'))),
-                         shape=len(ti6al4v.density_data_array))])
-# PointerType(create_type(np.float32)) PointerType(base_type=float32) PointerType(sp.Basic('float32'))
-
-Ti6Al4V_density_lookup = Functional_Parameter(
-    func='lookup',
-    args=[ti6al4v.density_temp_array_lookup, ti6al4v.density_data_array_lookup])
-
-Ti6Al4V_thermal_diffusivity1 = Solidus_Liquidus_Parameter(
-    sol=ti6al4v.thermal_diffusivity_solidus,
-    liq=ti6al4v.thermal_diffusivity_liquidus
-)
-
-Ti6Al4V_thermal_diffusivity = Functional_Parameter(
-    func='thermal_diffusivity',
-    args=[Ti6Al4V_heat_conductivity, Ti6Al4V_density_int, Ti6Al4V_specific_heat_capacity]
-)

src/materials/dataclasses.py

-from dataclasses import dataclass
-
-
-@dataclass
-class Solidus_Liquidus_Parameter:
-    sol: float
-    liq: float
-
-
-@dataclass
-class Functional_Parameter:
-    func: str
-    args: list

src/materials/elements.py

+from dataclasses import dataclass
+from src.materials.constants import Constants
+
+
+# Element Data Class
+@dataclass
+class Element:
+    atomic_number: float
+    atomic_mass: float
+    temperature_melt: float
+    temperature_boil: float
+    latent_heat_of_fusion: float
+    latent_heat_of_vaporization: float
+
+
+Ti = Element(
+    atomic_number=22,
+    atomic_mass=47.867*Constants.u,
+    temperature_melt=1941,
+    temperature_boil=3533,
+    latent_heat_of_fusion=18700,
+    latent_heat_of_vaporization=427000
+)
+
+
+Al = Element(
+    atomic_number=13,
+    atomic_mass=26.9815384*Constants.u,
+    temperature_melt=933.35,
+    temperature_boil=2743,
+    latent_heat_of_fusion=10700,
+    latent_heat_of_vaporization=284000
+)
+
+
+V = Element(
+    atomic_number=23,
+    atomic_mass=50.9415*Constants.u,
+    temperature_melt=2183,
+    temperature_boil=3680,
+    latent_heat_of_fusion=21500,
+    latent_heat_of_vaporization=444000
+)

src/materials/get_arguments.py

-# from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
-from get_parameters import get_parameters
-from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
-
-
-def get_arguments(T, temperatures, args):
-    """
-    Recursively get arguments for parameter functions.
-
-    :param T: Temperature value.
-    :param temperatures: Solidus and liquidus temperatures.
-    :param args: Arguments to be processed.
-    """
-    converted_args = []
-    for arg in args:
-        if isinstance(arg, (Solidus_Liquidus_Parameter, Functional_Parameter)):
-            converted_args.append(get_parameters(T, temperatures, arg))
-        else:
-            converted_args.append(arg)
-    return converted_args

src/materials/get_parameters.py

-import sympy as sp
-from src.materials.interpolators import interpolate, interpolate_property, density_lookup
-from src.materials.models import density_by_thermal_expansion, thermal_diffusivity_func
-# from src.materials.alloys.Ti6Al4V import Ti6Al4V_parameters
-from src.materials.dataclasses import Solidus_Liquidus_Parameter, Functional_Parameter
-
-
-def get_parameters(T, temperatures, parameter):
-    """
-    Resolve parameters based on temperature and parameter type.
-
-    :param T: Temperature value.
-    :param temperatures: Solidus and liquidus temperatures.
-    :param parameter: Parameter to be resolved.
-    """
-    if isinstance(parameter, Solidus_Liquidus_Parameter):
-        return sp.Piecewise(
-            (parameter.sol, T <= temperatures.sol),
-            (parameter.liq, T >= temperatures.liq),
-            (interpolate(T, temperatures.sol, temperatures.liq, parameter.sol, parameter.liq), True)
-        )
-    elif isinstance(parameter, Functional_Parameter):
-        resolved_args = [get_parameters(T, temperatures, arg) for arg in parameter.args]
-        if parameter.func == 'thermal_expansion':
-            return density_by_thermal_expansion(T, *resolved_args)  # *get_parameters(T, temperatures, parameter.args)
-        elif parameter.func == 'interpolation':
-            return interpolate_property(T, *resolved_args)
-        elif parameter.func == 'interpolation_in_array':  # TODO
-            T_base, T_incr, data_array, data_symbol = resolved_args
-            idx = sp.floor((T - T_base) / T_incr)
-            mod = (T - T_base) - idx * T_incr  # mod replaces (T - T_base) % T_incr
-            return sp.Piecewise(
-                (data_array[0], T < T_base),
-                (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
-                (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True)
-            )
-        elif parameter.func == 'lookup':
-            return density_lookup(T, *resolved_args)
-        elif parameter.func == 'thermal_diffusivity':
-            return thermal_diffusivity_func(*resolved_args)
-        else:
-            raise ValueError("No known format in the dict: parameter = ", parameter)
-    else:
-        return parameter  # assume everything which is not a dict is a float, int, etc.

src/materials/interpolators.py

 import numpy as np
 import sympy as sp
+from typing import Union, List
+# import pystencils as ps
+# from pystencils.types.quick import *
+from src.materials.property import Property
+from pystencils.sympyextensions import CastFunc
 
 
-def interpolate(x, xmin, xmax, ymin, ymax):
+def interpolate_equidistant(T: Union[float, sp.Symbol], T_base: float, T_incr: float,
+                            data_array: Union[np.ndarray, List[float]], label: str) -> Union[float, Property]:  # TODO
     """
-    Linear interpolation between two points: y = y_1 + ((y2 - y1) / (x2 - x1)) * (x - x1)
+    Interpolates values for equidistant temperature data points.
+
+    :param label:
+    :param T: Temperature value, either symbolic (sp.Expr) or numeric (float).
+    :param T_base: Base temperature value.
+    :param T_incr: Increment between successive temperature data points.
+    :param data_array: Array of data values corresponding to equidistant temperature points.
+    :return: Interpolated value at temperature T.
     """
-    return (ymin * (xmax - x) + ymax * (x - xmin)) / (xmax - xmin)
-
-
-def interpolate_property(T, base_temperature, base_values):  # TODO
+    if isinstance(T, sp.Symbol):
+        # Ensure all data_array elements are compatible with SymPy
+        # data_array = [sp.Float(float(x)) for x in data_array]
+        # data_symbol = sp.IndexedBase('data_symbol', shape=len(data_array))
+        # data_symbol = sp.IndexedBase(ps.TypedSymbol('density_data_symbol', Ptr(Fp(32))), shape=len(data_array))
+        data_array = tuple(sp.Float(float(x)) for x in data_array)  # Alternative
+        # data_symbol = ps.TypedSymbol(label, Arr(Fp(64)))  # Alternative
+        data_base = sp.IndexedBase(label, shape=len(data_array))  # Alternative
+        print(type(data_base))
+        print(type(data_base.label.name))
+        print(type(data_array))
+        # idx = sp.Symbol('idx', integer=True)
+        idx = sp.floor((T - T_base) / T_incr)
+        idx_int = CastFunc(idx, np.int64)
+        print(idx.is_Integer)
+        mod = (T - T_base) - idx * T_incr  # mod replaces (T - T_base) % T_incr
+        result = sp.Piecewise(
+            (data_array[0], T < T_base),
+            (data_array[-1], T >= T_base + (len(data_array) - 1) * T_incr),
+            # (((T_incr - mod) * data_symbol[idx] + (mod * data_symbol[idx + 1])) / T_incr, True))
+            (((T_incr - mod) * data_base[idx_int] + (mod * data_base[idx_int + 1])) / T_incr, True))  # Alternative
+        # return result.subs(data_symbol, sp.Array(data_array))
+        return Property(result, ((data_base, data_array),))  # [result, ps.Assignment(data_symbol, data_array), ps.Assignment(data_symbol, data_array)]  # Alternative resurn a property, 1st arg: result, 2nd: List[tuple]
+    else:
+        if T < T_base:
+            return data_array[0]
+        if T > T_base + (len(data_array) - 1) * T_incr:
+            return data_array[-1]
+        idx = int((T - T_base) / T_incr)
+        mod = (T - T_base) - idx * T_incr  # mod replaces (T - T_base) % T_incr
+        return ((T_incr - mod) * data_array[idx] + (mod * data_array[idx + 1])) / T_incr
+
+
+def interpolate_lookup(T: Union[float, sp.Symbol], base_temperature: Union[np.ndarray, List[float]],
+                       base_values: Union[np.ndarray, List[float]]) -> Union[float, sp.Expr]:  # TODO
     """
     Interpolates properties based on temperature for both symbolic and numeric cases.
 
@@ -17,98 +61,65 @@ def interpolate_property(T, base_temperature, base_values):  # TODO
     :param base_temperature: List or array of base temperatures.
     :param base_values: List or array of property values corresponding to base temperatures.
     """
-    # Ensure base_temperature and base_values are lists or arrays
-    if not isinstance(base_temperature, (list, np.ndarray)) or not isinstance(base_values, (list, np.ndarray)):
-        raise TypeError("base_temperature and base_values must be of type list or np.ndarray")
-    # Convert to lists if they are numpy arrays
-    if isinstance(base_temperature, np.ndarray):
-        base_temperature = base_temperature.tolist()
-    if isinstance(base_values, np.ndarray):
-        base_values = base_values.tolist()
-
-    # Debugging: Print lengths and types
-    # print(f"Length of base_temperature: {len(base_temperature)}")
-    # print(f"Length of base_values: {len(base_values)}")
-    # print(f"Type of base_temperature: {type(base_temperature)}")
-    # print(f"Type of base_values: {type(base_values)}")
+    # Ensure that base_temperature and base_values are numpy arrays or lists of floats
+    if isinstance(base_temperature, (np.ndarray, list)):
+        base_temperature = np.array(base_temperature, dtype=float)
+    if isinstance(base_values, (np.ndarray, list)):
+        base_values = np.array(base_values, dtype=float)
 
     # Check if base_temperature and base_values have the same length
     if len(base_temperature) != len(base_values):
         raise ValueError("base_temperature and base_values must have the same length")
     # assert len(base_temperature) == len(base_values), "base_temperature and base_values must have the same length"
 
-    if isinstance(T, sp.Expr):  # Symbolic expression
-        base_temperature = [sp.Float(t) for t in base_temperature]
-        base_values = [sp.Float(val) for val in base_values]
+    if isinstance(T, sp.Symbol):  # Symbolic expression
+        base_temperature = [sp.Float(float(t)) for t in base_temperature]
+        base_values = [sp.Float(float(val)) for val in base_values]
         conditions = [
             (base_values[0], T < base_temperature[0]),
             (base_values[-1], T >= base_temperature[-1])]
         for i in range(len(base_temperature) - 1):
             interp_expr = (base_values[i] +
-                           (base_values[i+1] - base_values[i]) / (base_temperature[i+1] - base_temperature[i]) *
+                           (base_values[i + 1] - base_values[i]) / (base_temperature[i + 1] - base_temperature[i]) *
                            (T - base_temperature[i]))
-            conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
+            if i + 1 == len(base_temperature) - 1:
+                conditions.append((interp_expr, sp.And(T >= base_temperature[i], T < base_temperature[i + 1])))
+            else:
+                conditions.append((interp_expr, True))
+        # return sp.Piecewise(*conditions)
         return sp.Piecewise(*conditions)
     else:  # Numeric expression
+        try:
+            T = float(T)  # Ensure T is a float for numeric case
+        except (TypeError, ValueError) as e:
+            raise ValueError(f"Invalid numeric input for T: {T}") from e
+
         return np.interp(T, base_temperature, base_values)
 
 
-# Example usage:
-Temp = sp.Symbol('T')
-density_temp_array1 = [200.0, 300.0, 400.0, 500.0, 600.0]
-density_data_array1 = [50.0, 40.0, 30.0, 20.0, 10.0]
-density_temp_array2 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
-density_data_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
+if __name__ == '__main__':
+    # Example usage:
+    Temp = sp.Symbol('Temp')
+    density_temp_array1 = np.array([200.0, 300.0, 400.0, 500.0, 600.0])
+    density_data_array1 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
+    density_temp_array2 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
+    density_data_array2 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
 
-# Symbolic interpolation
-symbolic_result = interpolate_property(Temp, density_temp_array2, density_data_array2)
-print("Symbolic interpolation result:", symbolic_result)
+    # Symbolic interpolation
+    symbolic_result = interpolate_lookup(Temp, density_temp_array1, density_data_array1)
+    print("Symbolic interpolation result with interpolate_lookup:", symbolic_result)
 
-# Numeric interpolation
-numeric_result = interpolate_property(450.0, density_temp_array2, density_data_array2)
-print("Numeric interpolation result:", numeric_result)
+    # Numeric interpolation
+    numeric_result = interpolate_lookup(450.0, density_temp_array1, density_data_array1)
+    print("Numeric interpolation result with interpolate_lookup:", numeric_result)
 
+    T_base = 200
+    T_incr = 100
 
-def density_lookup(T, temperature_array, density_array):
-    """
-    Look up density based on temperature using interpolation.
+    # Symbolic interpolation
+    symbolic_result = interpolate_equidistant(Temp, T_base, T_incr, density_data_array1, 'density')
+    print("Symbolic interpolation result with interpolate_equidistant:", symbolic_result)
 
-    :param T: Temperature value.
-    :param temperature_array: Array of temperatures.
-    :param density_array: Array of densities.
-    """
-    if len(temperature_array) != len(density_array):
-        raise ValueError("temperature_array and density_array must have the same length")
-    if isinstance(T, sp.Expr):  # Symbolic expression
-        conditions = [
-            (density_array[0], T < temperature_array[0]),
-            (density_array[-1], T >= temperature_array[-1])
-        ]
-        for i in range(len(temperature_array) - 1):
-            interp_expr = (density_array[i] +
-                           (density_array[i+1] - density_array[i]) / (temperature_array[i+1] - temperature_array[i]) *
-                           (T - temperature_array[i]))
-            conditions.append((interp_expr, sp.And(T >= temperature_array[i], T < temperature_array[i + 1])))
-        return sp.Piecewise(*conditions)
-    else:  # Numeric expression
-        return np.interp(T, temperature_array, density_array)
-
-
-# Example usage:
-Temp = sp.Symbol('T')
-density_temp_array3 = [200.0, 300.0, 400.0, 500.0, 600.0]
-density_data_array3 = [50.0, 40.0, 30.0, 20.0, 10.0]
-density_temp_array4 = np.array([200.0, 300.0, 400.0, 500.0, 600.0])
-density_data_array4 = np.array([50.0, 40.0, 30.0, 20.0, 10.0])
-density_temp_array5 = [1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0]
-density_data_array5 = [4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0]
-density_temp_array6 = np.array([1878.0, 1884.2, 1894.7, 1905.3, 1915.8, 1926.3, 1928.0])
-density_data_array6 = np.array([4236.3, 4234.9, 4233.6, 4232.3, 4230.9, 4229.7, 4227.0])
-
-# Symbolic interpolation
-symbolic_result = density_lookup(Temp, density_temp_array4, density_data_array4)
-print("Symbolic interpolation result:", symbolic_result)
-
-# Numeric interpolation
-numeric_result = density_lookup(450.0, density_temp_array4, density_data_array4)
-print("Numeric interpolation result:", numeric_result)
+    # Numeric interpolation
+    numeric_result = interpolate_equidistant(450.0, T_base, T_incr, density_data_array1, 'density')
+    print("Numeric interpolation result with interpolate_equidistant:", numeric_result)

src/materials/models.py

-def density_by_thermal_expansion(T, base_temperature, base_density, thermal_expansion_coefficient):
+from sympy import Expr, Symbol, simplify
+from typing import Union
+
+
+def density_by_thermal_expansion(
+        T: Union[float, Symbol],
+        base_temperature: float,
+        base_density: float,
+        thermal_expansion_coefficient: Union[float, Expr]) \
+        -> Union[float, Expr]:
     """
     Calculate density based on thermal expansion: rho_0 / (1 + tec * (T - T_0))
 
@@ -10,12 +19,19 @@ def density_by_thermal_expansion(T, base_temperature, base_density, thermal_expa
     return base_density / (1 + thermal_expansion_coefficient * (T - base_temperature)) ** 3
 
 
-def thermal_diffusivity_func(heat_conductivity, density, specific_heat_capacity):
+def thermal_diffusivity_by_heat_conductivity(
+        heat_conductivity: Union[float, Expr],
+        density: Union[float, Expr],
+        heat_capacity: Union[float, Expr]) \
+        -> Union[float, Expr]:
     """
     Calculate thermal diffusivity: k(T) / (rho(T) * c_p(T))
 
     :param heat_conductivity: Thermal conductivity.
     :param density: Density.
-    :param specific_heat_capacity: Specific heat capacity.
+    :param heat_capacity: Specific heat capacity.
     """
-    return heat_conductivity / (density * specific_heat_capacity)
+    result = heat_conductivity / (density * heat_capacity)
+    if isinstance(heat_conductivity, Expr) or isinstance(density, Expr) or isinstance(heat_capacity, Expr):
+        result = simplify(result)
+    return result

src/materials/property.py

+from sympy import Expr
+from dataclasses import dataclass
+
+
+@dataclass
+class Property:
+    expression: Expr | float
+    assignments: tuple[tuple, ...] = tuple