Skip to content
Snippets Groups Projects
Commit d41865b8 authored by Rahil Doshi's avatar Rahil Doshi
Browse files

Simplify energy_density to a single basic model

parent eeae9c50
Branches
No related tags found
No related merge requests found
Pipeline #77211 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
src/pymatlib/core/models.py
+ 1
− 238
View file @ d41865b8
...@@ -157,63 +157,7 @@ def thermal_diffusivity_by_heat_conductivity( ...@@ -157,63 +157,7 @@ def thermal_diffusivity_by_heat_conductivity(
raise ValueError("Division by zero encountered in thermal diffusivity calculation") raise ValueError("Division by zero encountered in thermal diffusivity calculation")
def energy_density_standard( def energy_density(
temperature: Union[float, sp.Symbol],
density: Union[float, MaterialProperty],
heat_capacity: Union[float, MaterialProperty],
latent_heat: Union[float, MaterialProperty]) \
-> MaterialProperty:
# Input validation to check for incompatible data types
# if isinstance(temperature, float) and temperature < ABSOLUTE_ZERO:
# raise ValueError(f"Temperature cannot be below absolute zero ({ABSOLUTE_ZERO}K)")
'''if isinstance(density, float) and density <= 0:
raise ValueError(f"Density must be positive, got {density}")
if isinstance(heat_capacity, float) and heat_capacity <= 0:
raise ValueError(f"Heat capacity must be positive, got {heat_capacity}")'''
# _validate_positive(density, heat_capacity, names=['Density', 'Heat capacity'])
'''if isinstance(latent_heat, float) and latent_heat < 0:
raise ValueError(f"Latent heat cannot be negative (should be zero or positive), got {latent_heat}")'''
# _validate_non_negative(latent_heat, names=['Latent heat'])
'''sub_assignments = [
assignment for prop in [density, heat_capacity, latent_heat]
if isinstance(prop, MaterialProperty)
for assignment in prop.assignments
]'''
(density_expr, heat_capacity_expr, latent_heat_expr), sub_assignments \
= _prepare_material_expressions(density, heat_capacity, latent_heat)
density_expr = density.expr if isinstance(density, MaterialProperty) else sympy_wrapper(density)
heat_capacity_expr = heat_capacity.expr if isinstance(heat_capacity, MaterialProperty) else sympy_wrapper(heat_capacity)
latent_heat_expr = latent_heat.expr if isinstance(latent_heat, MaterialProperty) else sympy_wrapper(latent_heat)
energy_density_expr = density_expr * (temperature * heat_capacity_expr + latent_heat_expr)
return MaterialProperty(energy_density_expr, sub_assignments)
# Just FYI: ps.Assignment(T.center, (s.h / s.density_mat - s.latent_heat_mat) / s.heat_capacity_mat)
# For backward compatibility
energy_density = energy_density_standard
def energy_density_enthalpy_based(
density: Union[float, MaterialProperty],
specific_enthalpy: Union[float, MaterialProperty],
latent_heat: Union[float, MaterialProperty]) \
-> MaterialProperty:
# _validate_positive(density, specific_enthalpy, names=["Density", "Specific enthalpy"])
# _validate_non_negative(latent_heat, names=["Latent heat"])
(density_expr, specific_enthalpy_expr, latent_heat_expr), sub_assignments \
= _prepare_material_expressions(density, specific_enthalpy, latent_heat)
energy_density_expr = density_expr * (specific_enthalpy_expr + latent_heat_expr)
return MaterialProperty(energy_density_expr, sub_assignments)
def energy_density_total_enthalpy(
density: Union[float, MaterialProperty], density: Union[float, MaterialProperty],
specific_enthalpy: Union[float, MaterialProperty]) \ specific_enthalpy: Union[float, MaterialProperty]) \
-> MaterialProperty: -> MaterialProperty:
...@@ -227,184 +171,3 @@ def energy_density_total_enthalpy( ...@@ -227,184 +171,3 @@ def energy_density_total_enthalpy(
#print(specific_enthalpy_expr) #print(specific_enthalpy_expr)
energy_density_expr = density_expr * specific_enthalpy_expr energy_density_expr = density_expr * specific_enthalpy_expr
return MaterialProperty(energy_density_expr, sub_assignments) return MaterialProperty(energy_density_expr, sub_assignments)
'''
def compute_energy_density_MaterialProperty(
temperature: Union[float, sp.Expr],
_density: Union[float, MaterialProperty],
heat_capacity: Union[float, MaterialProperty],
latent_heat: Union[float, MaterialProperty]) -> Tuple[MaterialProperty, float, float, float]:
"""
Compute energy density for MaterialProperty objects dynamically and return as a MaterialProperty.
"""
properties = {
"density": _density,
"heat_capacity": heat_capacity,
"latent_heat": latent_heat
}
evaluated_properties = {}
assignments = []
for name, prop in properties.items():
print(f"Processing {name}: {prop}, type: {type(prop)}")
if isinstance(prop, MaterialProperty):
print(f"{name} is a MaterialProperty")
if len(prop.assignments) == 0:
if isinstance(prop.expr, sp.Float):
evaluated_properties[name] = float(prop.expr)
elif isinstance(prop.expr, (sp.Piecewise, sp.Add, sp.Mul)):
free_symbols = prop.expr.free_symbols
if free_symbols:
substitutions = {symbol: temperature for symbol in free_symbols}
evaluated_properties[name] = prop.expr.subs(substitutions).evalf() if isinstance(temperature, float) else prop.expr.subs(substitutions)
else:
evaluated_properties[name] = prop.expr.evalf() if isinstance(temperature, float) else prop.expr
print(f"Evaluated {name}: {evaluated_properties[name]}")
elif len(prop.assignments) >= 1:
expr = prop.expr
subs_dict = {}
for assignment in prop.assignments:
ass_lhs = assignment.lhs
ass_rhs = assignment.rhs
print(f"{name} assignment - LHS: {ass_lhs}, RHS: {ass_rhs}")
if isinstance(ass_rhs, (sp.Piecewise, sp.Add)):
subs_dict_rhs = {symbol: temperature for symbol in ass_rhs.free_symbols}
ass_rhs_val = ass_rhs.subs(subs_dict_rhs).evalf() if isinstance(temperature, float) else ass_rhs.subs(subs_dict_rhs)
print(f"ass_rhs_val: {ass_rhs_val}")
ass_rhs_val = int(ass_rhs_val) if assignment.lhs_type == "int" and isinstance(temperature, float) else ass_rhs_val
subs_dict_lhs = {symbol: ass_rhs_val for symbol in ass_lhs.free_symbols}
subs_dict.update(subs_dict_rhs)
subs_dict.update(subs_dict_lhs)
print(f"updated subs_dict: {subs_dict}")
elif isinstance(ass_rhs, Tuple):
tuple_values = [
rhs_part.subs({symbol: temperature for symbol in rhs_part.free_symbols}).evalf()
if isinstance(temperature, float)
else rhs_part.subs({symbol: temperature for symbol in rhs_part.free_symbols})
for rhs_part in ass_rhs
]
subs_dict[ass_lhs] = ass_rhs
print(f"Updated subs_dict with tuple: {subs_dict}")
else:
raise TypeError(f"Unknown assignment type: {type(ass_rhs)}")
print(f"Substitution dictionary for {name}: {subs_dict}")
evaluated_properties[name] = expr.subs(subs_dict).evalf() if isinstance(temperature, float) else expr.subs(subs_dict)
print(f"Final evaluated {name}: {evaluated_properties[name]}")
assignments.extend(prop.assignments)
else:
raise NotImplementedError(f"Unsupported number of assignments for {name}: {len(prop.assignments)}")
else:
# Wrap the non-MaterialProperty inputs directly
evaluated_properties[name] = wrapper(prop)
# Compute the final energy density expression
energy_density_expr = evaluated_properties["density"] * (temperature * evaluated_properties["heat_capacity"] + evaluated_properties["latent_heat"])
# Combine all assignments into a MaterialProperty
energy_density_property = MaterialProperty(
expr=energy_density_expr,
assignments=assignments
)
print(f"Returning MaterialProperty with expr: {energy_density_property.expr} "
f"and assignments: {energy_density_property.assignments}")
return (
energy_density_property,
evaluated_properties["density"],
evaluated_properties["heat_capacity"],
evaluated_properties["latent_heat"],
)
def compute_energy_density_array(
density: Union[float, MaterialProperty],
heat_capacity: Union[float, MaterialProperty],
latent_heat: Union[float, MaterialProperty],
file_path: Optional[str] = None,
temperature_array: Optional[np.ndarray] = None) \
-> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Compute energy density for temperatures from a file or a NumPy array.
Args:
density (float or MaterialProperty): Density of the material.
heat_capacity (float or MaterialProperty): Specific heat capacity.
latent_heat (float or MaterialProperty): Latent heat of the material.
file_path (Optional[str]): Path to the temperature data file (optional).
temperature_array (Optional[np.ndarray]): Array of temperature values (optional).
Returns:
Tuple[np.ndarray, np.ndarray]:
- Array of temperatures.
- Array of corresponding energy densities.
Raises:
ValueError: If both `file_path` and `temperature_array` are provided, or neither is provided.
"""
print(f"Inside function\nfile_path: {file_path}, temperature_array: {temperature_array}")
# Validate input
if file_path is None and temperature_array is not None:
if isinstance(temperature_array, np.ndarray):
if temperature_array.size == 0:
raise ValueError(f"Temperature array must not be empty")
else:
raise TypeError(f"Temperature array must be a numpy array, got {type(temperature_array)}")
if file_path is not None and temperature_array is not None:
raise ValueError("Provide only one of `file_path` or `temperature_array`, not both.")
if file_path is None and temperature_array is None:
raise ValueError("You must provide either `file_path` or `temperature_array`.")
_temperatures = []
_densities = []
_heat_capacities = []
_latent_heats = []
_energy_densities = []
# Handle temperatures from file
if isinstance(file_path, str) and file_path:
try:
with open(file_path, 'r') as file:
for line in file:
# Skip non-data lines
if not line.strip() or not line[0].isdigit():
continue
# Extract the temperature from the first column
temperature = float(line.split()[0])
_temperatures.append(temperature)
except FileNotFoundError:
raise FileNotFoundError(f"The file '{file_path}' does not exist.")
except Exception as e:
raise ValueError(f"An error occurred while processing the file: {e}")
# Handle temperatures from NumPy array
if temperature_array is None:
raise ValueError("temperature_array is None")
if temperature_array is not None:
_temperatures = temperature_array.tolist()
print(f"_temperatures: {_temperatures}")
# Compute energy densities
for temperature in _temperatures:
print(f"temperature: {temperature}")
energy_density, density_val, specific_heat_val, latent_heat_val = compute_energy_density_MaterialProperty(
temperature=temperature,
_density=density,
heat_capacity=heat_capacity,
latent_heat=latent_heat,
)
_densities.append(density_val)
_heat_capacities.append(specific_heat_val)
_latent_heats.append(latent_heat_val)
_energy_densities.append(float(energy_density.expr))
# Convert lists to NumPy arrays
return (np.array(_temperatures),
np.array(_densities),
np.array(_heat_capacities),
np.array(_latent_heats),
np.array(_energy_densities))
'''
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment