diff --git a/apps/tutorials/codegen/HeatEquationKernel.py b/apps/tutorials/codegen/HeatEquationKernel.py
index 9c418fc977b34f24344619c04e40d565a2210980..a9c5e472ea25da75888b4bfa038ec175ffba6983 100644
--- a/apps/tutorials/codegen/HeatEquationKernel.py
+++ b/apps/tutorials/codegen/HeatEquationKernel.py
@@ -1,252 +1,11 @@
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)
diff --git a/src/materials/alloys/Alloy.py b/src/materials/alloys/Alloy.py
new file mode 100644
index 0000000000000000000000000000000000000000..27975f52fa1063e0a118b7c2a5990a57f5ae8d99
--- /dev/null
+++ b/src/materials/alloys/Alloy.py
@@ -0,0 +1,54 @@
+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)
diff --git a/src/materials/alloys/Ti6Al4V.py b/src/materials/alloys/Ti6Al4V.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b8ccd2ccfda5d31aecab5f6f8afd46deabb441d
--- /dev/null
+++ b/src/materials/alloys/Ti6Al4V.py
@@ -0,0 +1,65 @@
+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))
diff --git a/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py b/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py
deleted file mode 100644
index ff5897e704531a1e1244bd0c9053ccfa27b06100..0000000000000000000000000000000000000000
--- a/src/materials/alloys/Ti6Al4V/Ti6Al4V_parameters.py
+++ /dev/null
@@ -1,103 +0,0 @@
-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]
-)
diff --git a/src/materials/alloys/Ti6Al4V/__init__.py b/src/materials/alloys/Ti6Al4V/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/materials/dataclasses.py b/src/materials/dataclasses.py
deleted file mode 100644
index a67834b2d70628044a2da0986dbed7731c23b4a2..0000000000000000000000000000000000000000
--- a/src/materials/dataclasses.py
+++ /dev/null
@@ -1,13 +0,0 @@
-from dataclasses import dataclass
-
-
-@dataclass
-class Solidus_Liquidus_Parameter:
- sol: float
- liq: float
-
-
-@dataclass
-class Functional_Parameter:
- func: str
- args: list
diff --git a/src/materials/elements.py b/src/materials/elements.py
new file mode 100644
index 0000000000000000000000000000000000000000..049fc1756d98550a2f56042807f88124cc66608b
--- /dev/null
+++ b/src/materials/elements.py
@@ -0,0 +1,43 @@
+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
+)
diff --git a/src/materials/elements/Ti.py b/src/materials/elements/Ti.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/materials/elements/__init__.py b/src/materials/elements/__init__.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/src/materials/get_arguments.py b/src/materials/get_arguments.py
deleted file mode 100644
index d6a3bc2e0aca0c3be9ba58ac4b9d87ef2dfcc602..0000000000000000000000000000000000000000
--- a/src/materials/get_arguments.py
+++ /dev/null
@@ -1,20 +0,0 @@
-# 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
diff --git a/src/materials/get_parameters.py b/src/materials/get_parameters.py
deleted file mode 100644
index 311ac0eda1cf9470340166e322c26aa3be0a3868..0000000000000000000000000000000000000000
--- a/src/materials/get_parameters.py
+++ /dev/null
@@ -1,44 +0,0 @@
-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.
diff --git a/src/materials/interpolators.py b/src/materials/interpolators.py
index 242254280a2e6bf712862ac077744f67c22fda12..6b32bf27768b537930b26535f92f6c28cb9c5808 100644
--- a/src/materials/interpolators.py
+++ b/src/materials/interpolators.py
@@ -1,15 +1,59 @@
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)
diff --git a/src/materials/models.py b/src/materials/models.py
index cd098fec05fdc4899f95e48450874ce3386d1a8f..a4a7372eae51d65e9d8454623148c032d962f480 100644
--- a/src/materials/models.py
+++ b/src/materials/models.py
@@ -1,4 +1,13 @@
-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
diff --git a/src/materials/property.py b/src/materials/property.py
new file mode 100644
index 0000000000000000000000000000000000000000..aacd51db2d94a68084ecae3a8545c36a2ff1b6a4
--- /dev/null
+++ b/src/materials/property.py
@@ -0,0 +1,8 @@
+from sympy import Expr
+from dataclasses import dataclass
+
+
+@dataclass
+class Property:
+ expression: Expr | float
+ assignments: tuple[tuple, ...] = tuple