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

Merge branch 'development' into release/v0.1.0

parents 6e1e556e e422044f
Branches release/v0.1.0
No related merge requests found
Pipeline #77633 passed with stage
in 22 seconds
import re
import numpy as np
from typing import List, Any
from pystencils.types import PsCustomType
from pystencilssfg import SfgComposer
from pystencilssfg.composer.custom import CustomGenerator
......@@ -7,55 +8,9 @@ from pymatlib.core.interpolators import prepare_interpolation_arrays
class InterpolationArrayContainer(CustomGenerator):
"""Container for x-y interpolation arrays and methods.
This class stores x and y arrays and generates C++ code
for efficient conversion to compute y for a given x. It supports both
binary search interpolation (O(log n)) and double lookup interpolation (O(1))
with automatic method selection based on data characteristics.
Attributes:
name (str): Name for the generated C++ class.
x_array (np.ndarray): Array of x values (must be monotonically increasing).
y_array (np.ndarray): Array of y values corresponding to x_array.
method (str): Interpolation method selected ("binary_search" or "double_lookup").
x_bs (np.ndarray): x array prepared for binary search.
y_bs (np.ndarray): y array prepared for binary search.
has_double_lookup (bool): Whether double lookup interpolation is available.
If has_double_lookup is True, the following attributes are also available:
x_eq (np.ndarray): Equidistant x array for double lookup.
y_neq (np.ndarray): Non-equidistant y array for double lookup.
y_eq (np.ndarray): Equidistant y array for double lookup.
inv_delta_y_eq (float): Inverse of the y step size for double lookup.
idx_map (np.ndarray): Index mapping array for double lookup.
Examples:
>>> import numpy as np
>>> from pystencilssfg import SfgComposer
>>> from pymatlib.core.codegen.interpolation_array_container import InterpolationArrayContainer
>>>
>>> # Create temperature and energy arrays
>>> T = np.array([300, 600, 900, 1200], dtype=np.float64)
>>> E = np.array([1e9, 2e9, 3e9, 4e9], dtype=np.float64)
>>>
>>> # Create and generate the container
>>> with SfgComposer() as sfg:
>>> container = InterpolationArrayContainer("MyMaterial", T, E)
>>> sfg.generate(container)
"""
"""Container for x-y interpolation arrays and methods."""
def __init__(self, name: str, x_array: np.ndarray, y_array: np.ndarray):
"""Initialize the interpolation container.
Args:
name (str): Name for the generated C++ class.
x_array (np.ndarray): Array of x values.
Must be monotonically increasing.
y_array (np.ndarray): Array of y values
corresponding to x_array.
Raises:
ValueError: If arrays are empty, have different lengths, or are not monotonic.
TypeError: If name is not a string or arrays are not numpy arrays.
"""
"""Initialize the interpolation container."""
super().__init__()
# Validate inputs
......@@ -66,7 +21,7 @@ class InterpolationArrayContainer(CustomGenerator):
raise ValueError(f"'{name}' is not a valid C++ class name")
if not isinstance(x_array, np.ndarray) or not isinstance(y_array, np.ndarray):
raise TypeError("Temperature and energy arrays must be numpy arrays")
raise TypeError("x_array and y_array must be numpy arrays")
self.name = name
self.x_array = x_array
......@@ -95,12 +50,7 @@ class InterpolationArrayContainer(CustomGenerator):
raise ValueError(f"Failed to prepare interpolation arrays: {e}") from e
def _generate_binary_search(self, sfg: SfgComposer):
"""Generate code for binary search interpolation.
Args:
sfg (SfgComposer): Source file generator composer.
Returns:
list: List of public members for the C++ class.
"""
"""Generate code for binary search interpolation."""
x_bs_arr_values = ", ".join(str(v) for v in self.x_bs)
y_bs_arr_values = ", ".join(str(v) for v in self.y_bs)
......@@ -117,13 +67,8 @@ class InterpolationArrayContainer(CustomGenerator):
)
]
def _generate_double_lookup(self, sfg: SfgComposer):
"""Generate code for double lookup interpolation.
Args:
sfg (SfgComposer): Source file generator composer.
Returns:
list: List of public members for the C++ class.
"""
def _generate_double_lookup(self, sfg: SfgComposer) -> List[Any]:
"""Generate code for double lookup interpolation."""
if not self.has_double_lookup:
return []
......@@ -148,26 +93,19 @@ class InterpolationArrayContainer(CustomGenerator):
)
]
def generate(self, sfg: SfgComposer):
"""Generate C++ code for the interpolation container.
This method generates a C++ class with the necessary arrays and methods
for temperature-energy interpolation.
Args:
sfg (SfgComposer): Source file generator composer.
"""
def generate(self, sfg: SfgComposer) -> None:
"""Generate C++ code for the interpolation container."""
sfg.include("<array>")
sfg.include("pymatlib_interpolators/interpolate_binary_search_cpp.h")
public_members = self._generate_binary_search(sfg)
# Add double lookup if available
# Add interpolate method that uses recommended approach
y_target = sfg.var("y_target", "double")
if self.has_double_lookup:
sfg.include("pymatlib_interpolators/interpolate_double_lookup_cpp.h")
public_members.extend(self._generate_double_lookup(sfg))
# Add interpolate method that uses recommended approach
y_target = sfg.var("y_target", "double")
if self.has_double_lookup:
public_members.append(
sfg.method("interpolate", returns=PsCustomType("[[nodiscard]] double"), inline=True, const=True)(
sfg.expr("return interpolate_double_lookup_cpp({}, *this);", y_target)
......
......@@ -33,12 +33,13 @@ double interpolate_binary_search_cpp(
}
}
// Linear interpolation
const double x1 = arrs.y_bs[right];
const double x2 = arrs.y_bs[left];
const double y1 = arrs.x_bs[right];
const double y2 = arrs.x_bs[left];
// Get interpolation points
const double x1 = arrs.x_bs[left];
const double x2 = arrs.x_bs[right];
const double y1 = arrs.y_bs[left];
const double y2 = arrs.y_bs[right];
const double slope = (y2 - y1) / (x2 - x1);
return y1 + slope * (y_target - x1);
// Linear interpolation
const double inv_slope = (x2 - x1) / (y2 - y1);
return x1 + inv_slope * (y_target - y1);
}
......@@ -26,12 +26,12 @@ double interpolate_double_lookup_cpp(
idx_y_neq += arrs.y_neq[idx_y_neq + 1] < y_target;
// Get interpolation points
const double y1 = arrs.y_neq[idx_y_neq];
const double y2 = arrs.y_neq[idx_y_neq + 1];
const double x1 = arrs.x_eq[idx_y_neq];
const double x2 = arrs.x_eq[idx_y_neq + 1];
const double y1 = arrs.y_neq[idx_y_neq];
const double y2 = arrs.y_neq[idx_y_neq + 1];
// Linear interpolation
const double slope = (x2 - x1) / (y2 - y1);
return x1 + slope * (y_target - y1);
const double inv_slope = (x2 - x1) / (y2 - y1);
return x1 + inv_slope * (y_target - y1);
}
from dataclasses import dataclass
@dataclass
class Constants:
"""
A dataclass to store fundamental constants.
Attributes:
temperature_room (float): Room temperature in Kelvin.
N_a (float): Avogadro's number, the number of constituent particles (usually atoms or molecules) in one mole of a given substance.
amu (float): Atomic mass unit in kilograms.
e (float): Elementary charge, the electric charge carried by a single proton or the magnitude of the electric charge carried by a single electron.
speed_of_light (float): Speed of light in vacuum in meters per second.
TEMPERATURE_ROOM (float): Room temperature in Kelvin.
N_A (float): Avogadro's number, the number of constituent particles (usually atoms or molecules) in one mole of a given substance.
AMU (float): Atomic mass unit in kilograms.
E (float): Elementary charge, the electric charge carried by a single proton or the magnitude of the electric charge carried by a single electron.
SPEED_OF_LIGHT (float): Speed of light in vacuum in meters per second.
"""
temperature_room: float = 298.15 # Room temperature in Kelvin
N_a: float = 6.022141e23 # Avogadro's number, /mol
amu: float = 1.660538e-27 # Atomic mass unit, kg
e: float = 1.60217657e-19 # Elementary charge, C
speed_of_light: float = 0.299792458e9 # Speed of light, m/s
TEMPERATURE_ROOM: float = 298.15 # Room temperature in Kelvin
N_A: float = 6.022141e23 # Avogadro's number, /mol
AMU: float = 1.660538e-27 # Atomic mass unit, kg
E: float = 1.60217657e-19 # Elementary charge, C
SPEED_OF_LIGHT: float = 0.299792458e9 # Speed of light, m/s
......@@ -5,130 +5,130 @@ from pymatlib.core.elements import ChemicalElement
# RSC: Royal Society of Chemistry
# CRC: CRC Handbook of Chemistry and Physics
C = ChemicalElement(
CARBON = ChemicalElement(
name="Carbon",
atomic_number=6,
atomic_mass=12.0107 * Constants.amu,
atomic_mass=12.0107 * Constants.AMU,
temperature_melt=3915, # Melting temperature = 3915 K
temperature_boil=4300, # Boiling temperature = 4300 K
latent_heat_of_fusion=117000, # Latent heat of fusion = 117 kJ/mol
latent_heat_of_vaporization=355000 # Latent heat of vaporization = 355 kJ/mol
)
N = ChemicalElement(
NITROGEN = ChemicalElement(
name="Nitrogen",
atomic_number=7,
atomic_mass=14.0067 * Constants.amu,
atomic_mass=14.0067 * Constants.AMU,
temperature_melt=63.15, # Melting temperature = 63.15 K
temperature_boil=77.36, # Boiling temperature = 77.36 K
latent_heat_of_fusion=720, # Latent heat of fusion = 0.72 kJ/mol
latent_heat_of_vaporization=5570 # Latent heat of vaporization = 5.57 kJ/mol
)
Al = ChemicalElement(
ALUMINIUM = ChemicalElement(
name="Aluminium",
atomic_number=13, # Atomic number = 13 / Source: Periodic Table
atomic_mass=26.9815384 * Constants.amu, # Atomic mass = 26.9815384 amu / Source: NIST
atomic_mass=26.9815384 * Constants.AMU, # Atomic mass = 26.9815384 amu / Source: NIST
temperature_melt=933.35, # Melting temperature = 933.35 K / Source: RSC
temperature_boil=2743, # Boiling temperature = 2743 K / Source: RSC
latent_heat_of_fusion=10700, # Latent heat of fusion = 10700 J/kg / Source: CRC
latent_heat_of_vaporization=284000 # Latent heat of vaporization = 284000 J/kg / Source: CRC
)
Si = ChemicalElement(
SILICON = ChemicalElement(
name="Silicon",
atomic_number=14,
atomic_mass=28.0855 * Constants.amu,
atomic_mass=28.0855 * Constants.AMU,
temperature_melt=1687, # Melting temperature = 1687 K
temperature_boil=3538, # Boiling temperature = 3538 K
latent_heat_of_fusion=50200, # Latent heat of fusion = 50.2 kJ/mol
latent_heat_of_vaporization=359000 # Latent heat of vaporization = 359 kJ/mol
)
P = ChemicalElement(
PHOSPHORUS = ChemicalElement(
name="Phosphorus",
atomic_number=15,
atomic_mass=30.973762 * Constants.amu,
atomic_mass=30.973762 * Constants.AMU,
temperature_melt=317.3, # Melting temperature = 317.3 K
temperature_boil=553.7, # Boiling temperature = 553.7 K
latent_heat_of_fusion=2510, # Latent heat of fusion = 2.51 kJ/mol
latent_heat_of_vaporization=12400 # Latent heat of vaporization = 12.4 kJ/mol
)
S = ChemicalElement(
SULFUR = ChemicalElement(
name="Sulfur",
atomic_number=16,
atomic_mass=32.065 * Constants.amu,
atomic_mass=32.065 * Constants.AMU,
temperature_melt=388.36, # Melting temperature = 388.36 K
temperature_boil=717.8, # Boiling temperature = 717.8 K
latent_heat_of_fusion=1730, # Latent heat of fusion = 1.73 kJ/mol
latent_heat_of_vaporization=9800 # Latent heat of vaporization = 9.8 kJ/mol
)
Ti = ChemicalElement(
TITANIUM = ChemicalElement(
name="Titanium",
atomic_number=22, # Atomic number = 22 / Source: Periodic Table
atomic_mass=47.867 * Constants.amu, # Atomic mass = 47.867 amu / Source: NIST
atomic_mass=47.867 * Constants.AMU, # Atomic mass = 47.867 amu / Source: NIST
temperature_melt=1941, # Melting temperature = 1941 K / Source: RSC
temperature_boil=3560, # Boiling temperature = 3560 K / Source: RSC
latent_heat_of_fusion=18700, # Latent heat of fusion = 18700 J/kg / Source: CRC
latent_heat_of_vaporization=427000 # Latent heat of vaporization = 427000 J/kg / Source: CRC
)
V = ChemicalElement(
VANADIUM = ChemicalElement(
name="Vanadium",
atomic_number=23, # Atomic number = 23 / Source: Periodic Table
atomic_mass=50.9415 * Constants.amu, # Atomic mass = 50.9415 amu / Source: NIST
atomic_mass=50.9415 * Constants.AMU, # Atomic mass = 50.9415 amu / Source: NIST
temperature_melt=2183, # Melting temperature = 2183 K / Source: RSC
temperature_boil=3680, # Boiling temperature = 3680 K / Source: RSC
latent_heat_of_fusion=21500, # Latent heat of fusion = 21500 J/kg / Source: CRC
latent_heat_of_vaporization=444000 # Latent heat of vaporization = 444000 J/kg / Source: CRC
)
Cr = ChemicalElement(
CHROMIUM = ChemicalElement(
name="Chromium",
atomic_number=24, # Atomic number = 24 / Source: Periodic Table
atomic_mass=51.9961 * Constants.amu, # Atomic mass = 51.9961 amu / Source: NIST
atomic_mass=51.9961 * Constants.AMU, # Atomic mass = 51.9961 amu / Source: NIST
temperature_melt=2180, # Melting temperature = 2180 K / Source: RSC
temperature_boil=2944, # Boiling temperature = 2944 K / Source: RSC
latent_heat_of_fusion=16500, # Latent heat of fusion = 16500 J/kg / Source: CRC
latent_heat_of_vaporization=344000 # Latent heat of vaporization = 344000 J/kg / Source: CRC
)
Mn = ChemicalElement(
MANGANESE = ChemicalElement(
name="Manganese",
atomic_number=25, # Atomic number = 25 / Source: Periodic Table
atomic_mass=54.938045 * Constants.amu, # Atomic mass = 54.938045 amu / Source: NIST
atomic_mass=54.938045 * Constants.AMU, # Atomic mass = 54.938045 amu / Source: NIST
temperature_melt=1519, # Melting temperature = 1519 K / Source: RSC
temperature_boil=2334, # Boiling temperature = 2334 K / Source: RSC
latent_heat_of_fusion=12500, # Latent heat of fusion = 12500 J/kg / Source: CRC
latent_heat_of_vaporization=220000 # Latent heat of vaporization = 220000 J/kg / Source: CRC
)
Fe = ChemicalElement(
IRON = ChemicalElement(
name="Iron",
atomic_number=26, # Atomic number = 26 / Source: Periodic Table
atomic_mass=55.845 * Constants.amu, # Atomic mass = 55.845 amu / Source: NIST
atomic_mass=55.845 * Constants.AMU, # Atomic mass = 55.845 amu / Source: NIST
temperature_melt=1809, # Melting temperature = 1809 K / Source: RSC
temperature_boil=3134, # Boiling temperature = 3134 K / Source: RSC
latent_heat_of_fusion=13800, # Latent heat of fusion = 13800 J/kg / Source: CRC
latent_heat_of_vaporization=340000 # Latent heat of vaporization = 340000 J/kg / Source: CRC
)
Ni = ChemicalElement(
NICKEL = ChemicalElement(
name="Nickel",
atomic_number=28, # Atomic number = 28 / Source: Periodic Table
atomic_mass=58.6934 * Constants.amu, # Atomic mass = 58.6934 amu / Source: NIST
atomic_mass=58.6934 * Constants.AMU, # Atomic mass = 58.6934 amu / Source: NIST
temperature_melt=1728, # Melting temperature = 1728 K / Source: RSC
temperature_boil=3186, # Boiling temperature = 3186 K / Source: RSC
latent_heat_of_fusion=17200, # Latent heat of fusion = 17200 J/kg / Source: CRC
latent_heat_of_vaporization=377000 # Latent heat of vaporization = 377000 J/kg / Source: CRC
)
Mo = ChemicalElement(
MOLYBDENUM = ChemicalElement(
name="Molybdenum",
atomic_number=42,
atomic_mass=95.96 * Constants.amu,
atomic_mass=95.96 * Constants.AMU,
temperature_melt=2896, # Melting temperature = 2896K (2623°C)
temperature_boil=4912, # Boiling temperature = 4912K (4639°C)
latent_heat_of_fusion=37480, # Latent heat of fusion = 37.48 kJ/mol
......@@ -138,17 +138,17 @@ Mo = ChemicalElement(
# This dictionary maps chemical symbols (strings) to their corresponding ChemicalElement instances,
# allowing the parser to convert composition keys from the YAML file (like 'Fe': 0.675) to actual ChemicalElement objects needed by the Alloy class.
element_map = {
'C': C,
'N': N,
'Al': Al,
'Si': Si,
'P': P,
'S': S,
'Ti': Ti,
'V': V,
'Cr': Cr,
'Mn': Mn,
'Fe': Fe,
'Ni': Ni,
'Mo': Mo,
'C': CARBON,
'N': NITROGEN,
'Al': ALUMINIUM,
'Si': SILICON,
'P': PHOSPHORUS,
'S': SULFUR,
'Ti': TITANIUM,
'V': VANADIUM,
'Cr': CHROMIUM,
'Mn': MANGANESE,
'Fe': IRON,
'Ni': NICKEL,
'Mo': MOLYBDENUM,
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment