Add docs

docs/explanation/design_philosophy.md

+# Design Philosophy of pymatlib
+
+This document explains the core design principles, architectural decisions, and the rationale behind pymatlib's structure and implementation.
+
+## Core Principles
+
+pymatlib is built upon several core principles:
+
+- **Modularity**: Clearly separated components for ease of maintenance, testing, and extensibility.
+- **Flexibility**: Allow users to define material properties in various intuitive ways.
+- **Performance**: Leverage symbolic computation and optimized C++ code generation for high-performance simulations.
+- **Transparency and Reproducibility**: Clearly document material property definitions and computations to ensure reproducibility.
+
+## Layered Architecture
+
+pymatlib follows a layered architecture to separate concerns clearly:
+
+1. User Interface Layer (YAML Configuration)
+
+- Provides a simple, human-readable format for defining alloys and their properties.
+- Allows users to specify properties using multiple intuitive methods (constants, interpolation points, file-based data, computed properties).
+- Ensures clarity, readability, and ease of use
+
+2. Symbolic Representation Layer (SymPy)
+
+- Uses symbolic mathematics (via SymPy) internally to represent material properties.
+- Enables symbolic manipulation, simplification, and validation of property definitions.
+- Facilitates automatic computation of derived properties.
+
+3. Interpolation and Computation Layer (Python)
+
+- Implements robust interpolation methods (`interpolate_equidistant`, `interpolate_lookup`) for evaluating temperature-dependent properties.
+- Automatically analyzes input arrays (`prepare_interpolation_arrays`) to determine the optimal interpolation method based on data characteristics.
+- Provides symbolic manipulation capabilities through SymPy for computed properties.
+
+4. Code Generation Layer (C++)
+
+- Uses InterpolationArrayContainer to generate optimized C++ code for efficient energy-to-temperature conversions.
+- Automatically selects the best interpolation method (Binary Search or Double Lookup) based on data analysis results from Python.
+- Integrates seamlessly with pystencils and waLBerla for high-performance simulations.
+
+## Separation Between Python and C++ Components
+
+pymatlib separates responsibilities clearly between Python and C++ components:
+
+| Component                      | Responsibility                                                  | Implementation Language |
+|--------------------------------|--------------------------------------------------|----------|
+| YAML Configuration             | User-friendly material property definitions     | YAML     |
+| Symbolic Computation & Validation| Parsing YAML files, symbolic computations and validations (sympy) | Python   |
+| Interpolation Analysis           | Determining optimal interpolation method based on data characteristics (`prepare_interpolation_arrays`)       | Python   |
+| Code Generation                  | Generating optimized interpolation code (`InterpolationArrayContainer`) for high-performance simulations       | C++      |
+| Simulation Execution             | Running generated kernels within simulation frameworks (waLBerla/pystencils)   | C++      |
+
+This separation ensures:
+
+- Flexibility in defining materials through easily editable YAML files
+- Powerful symbolic manipulation capabilities in Python
+- High-performance numerical computations in C++
+
+## Why YAML?
+
+YAML was chosen as the primary configuration format because:
+
+- It is human-readable and easy to edit manually
+- It naturally supports nested structures required by complex material definitions
+- It integrates smoothly with Python's ecosystem via libraries like PyYAML
+- It allows referencing previously defined variables within the file (e.g., `solidus_temperature`, `liquidus_temperature`), reducing redundancy.
+
+## Integration with pystencils and waLBerla
+
+pymatlib integrates closely with [pystencils](https://pycodegen.pages.i10git.cs.fau.de/pystencils/) and [waLBerla](https://walberla.net/) through the following workflow:
+
+1. **Symbolic Definition**: Material properties are defined symbolically in pymatlib using YAML configurations.
+2. **Assignment Conversion**: The assignment_converter function converts pymatlib's symbolic assignments into pystencils-compatible assignments.
+3. **Code Generation**: pystencils generates optimized kernels from these assignments for numerical simulations.
+4. **Simulation execution**: Generated kernels are executed within waLBerla frameworks for large-scale parallel simulations.
+
+This integration allows pymatlib to leverage:
+
+- Symbolic mathematics from sympy via pystencils
+- Optimized stencil-based numerical kernels generated by pystencils-sfg
+- High-performance parallel computing capabilities provided by waLBerla
+
+## Automatic Method Selection for Interpolation
+
+A key design decision in pymatlib is automatic method selection for interpolation between energy density and temperature:
+
+1. **Analysis Phase (Python)**: 
+
+- The function prepare_interpolation_arrays() analyzes input arrays to determine if they're equidistant or not.
+- Based on this analysis, it selects either Binary Search or Double Lookup as the preferred interpolation method.
+
+2. **Code Generation Phase (C++)**:
+
+- The InterpolationArrayContainer class generates optimized C++ code that includes both interpolation methods (interpolateBS, interpolateDL) if applicable.
+- A wrapper method (interpolate) is automatically generated to select the best available method at runtime without user intervention.
+
+This ensures optimal performance without burdening users with manual selection decisions.
+
+## Extensibility
+
+pymatlib is designed with extensibility in mind:
+
+- Users can easily define new material properties or extend existing ones through YAML files without changing core code.
+- New computational models or interpolation methods can be added at the Python layer without affecting existing functionality.
+- The modular design allows easy integration with additional simulation frameworks beyond pystencils or waLBerla if needed.
+
+## Robustness and Error Handling
+
+pymatlib includes comprehensive validation checks at every step:
+
+### YAML Parser Validation
+
+- Ensures consistent units (SI units recommended).
+- Checks monotonicity of temperature-energy arrays.
+- Validates dependencies among computed properties.
+
+### Interpolation Validation
+
+The interpolation functions validate:
+
+- Array lengths match
+- Arrays contain sufficient elements
+- Arrays are strictly monotonic
+- Energy density increases consistently with temperature
+
+If any validation fails, clear error messages guide users toward correcting their configurations.
+
+## Performance Optimization Philosophy
+
+Performance-critical numerical operations are implemented in optimized C++ code generated automatically from Python definitions. This approach provides:
+
+1. **Ease of use**: Users define materials symbolically or numerically without worrying about low-level optimization details.
+2. **Automatic Optimization**: pymatlib automatically selects optimal algorithms (Double Lookup vs Binary Search) based on data characteristics without user intervention.
+3. **High Performance**: Generated kernels provide near-native performance suitable for large-scale simulations while maintaining flexibility in material definitions.
+

docs/explanation/interpolation_methods.md

+# Interpolation Methods in pymatlib
+
+This document explains the interpolation techniques used in pymatlib for energy-temperature conversion,
+including their internal workings, automatic method selection, and performance considerations.
+
+### Overview
+
+Interpolation is essential in pymatlib for converting between energy density and temperature values, particularly when modeling materials whose properties vary significantly with temperature. pymatlib provides two interpolation methods:
+
+- Binary Search Interpolation
+- Double Lookup Interpolation
+
+Additionally, pymatlib offers an intelligent wrapper method that automatically selects the optimal interpolation method based on data characteristics.
+
+## Binary Search Interpolation (`interpolateBS`)
+
+Binary Search Interpolation is one of the primary methods used in pymatlib for obtaining temperatures from energy density values.
+
+### How It Works
+
+1. The algorithm starts with a sorted array of energy density values and their corresponding temperature values
+2. When given an energy density value to convert to temperature:
+- It performs a binary search to find the position of the target value
+- It divides the search interval in half repeatedly until finding the closest match
+- Once the position is found, it performs linear interpolation between the two closest points
+
+### Implementation Details
+
+The binary search implementation handles edge cases gracefully:
+- If the target energy is below the minimum value, it returns the lowest temperature
+- If the target energy is above the maximum value, it returns the highest temperature
+- For values within the range, it calculates the interpolated temperature using:
+```text
+T = T[i] + (T[i+1] - T[i]) * (E_target - E[i]) / (E[i+1] - E[i])
+```
+
+As shown in the test code, this method is robust and handles a wide range of inputs, including extreme values and edge cases.
+
+### Characteristics
+- **Time Complexity**: O(log n) where n is the number of data points
+- **Advantages**: Works reliably with any monotonically increasing data
+- **Implementation**: Available through the `interpolateBS()` method in generated C++ code
+
+## Double Lookup Interpolation (`interpolateDL`)
+
+Double Lookup Interpolation is an optimized approach designed specifically for uniform or near-uniform data distributions.
+
+### How It Works
+
+1. During initialization (Python)
+- The algorithm checks if the temperature array is equidistant using `check_equidistant()`
+- It creates a uniform grid of energy values (`E_eq`) with carefully calculated spacing:
+```python
+delta_min = np.min(np.abs(np.diff(E_neq)))
+delta_E_eq = max(np.floor(delta_min * 0.95), 1.)
+E_eq = np.arange(E_neq[0], E_neq[-1] + delta_E_eq, delta_E_eq)
+```
+- It pre-computes an index mapping between the original energy array and the uniform grid:
+```python
+idx_map = np.searchsorted(E_neq, E_eq, side='right') - 1
+idx_map = np.clip(idx_map, 0, len(E_neq) - 2)
+```
+- It calculates the inverse of the energy grid spacing for fast index calculation:
+
+```python
+inv_delta_E_eq = 1.0 / (E_eq[1] - E_eq[0])
+```
+
+2. During lookup (C++):
+- It calculates the array index directly (O(1)) using the pre-computed inverse delta:
+```python
+idx = floor((E_target - E_eq[0]) * inv_delta_E_eq)
+```
+- It retrieves the pre-computed mapping to find the correct segment in the original array
+- It performs a simple linear interpolation between the two closest points
+
+### Characteristics
+- **Time Complexity**: O(1) - constant time regardless of array size using pre-computed mappings
+- **Best For**: Uniform temperature distributions
+- **Advantages**: Significantly faster (typically 2-4x) than binary search for uniform data
+- **Implementation**: Available through the `interpolateDL()` method in generated C++ code
+
+### Performance Optimization
+
+The double lookup method achieves O(1) complexity through several optimizations:
+- Pre-computing the index mapping eliminates the need for binary search
+- Using the inverse of the delta for multiplication instead of division
+- Constraining the index mapping to valid bounds during initialization
+
+## Automatic Method Selection Process (`interpolate`)
+
+```c++
+constexpr double energy_density = 16950000000.0;
+constexpr SimpleSteel alloy;
+const double temperature = alloy.interpolate(energy_density);
+```
+
+How Automatic Selection Works Internally
+
+The automatic method selection in pymatlib involves two key components working together:
+
+1. **Analysis Phase (Python)**: The prepare_interpolation_arrays() function analyzes the input temperature and energy density arrays to determine the most appropriate interpolation method:
+- It validates that both arrays are monotonic and energy increases with temperature
+- It checks if the temperature array is equidistant using check_equidistant()
+- Based on this analysis, it selects either "binary_search" or "double_lookup" as the method
+- It returns a dictionary containing the processed arrays and the selected method
+
+```python
+if is_equidistant:
+    E_eq, inv_delta_E_eq = E_eq_from_E_neq(E_bs)
+    idx_mapping = create_idx_mapping(E_bs, E_eq)
+    # Set method to "double_lookup"
+else:
+    # Set method to "binary_search"
+```
+
+2. **Code Generation Phase (C++)**: The InterpolationArrayContainer class uses this information to generate optimized C++ code:
+- It always includes the binary search method (interpolateBS)
+- If the temperature array is equidistant, it also includes the double lookup method (interpolateDL)
+- It adds an intelligent wrapper method interpolate() that automatically calls the preferred method:
+
+```c++
+// If double lookup is available, use it
+if (self.has_double_lookup):
+    interpolate() { return interpolate_double_lookup_cpp(E_target, *this); }
+// Otherwise, fall back to binary search
+else:
+    interpolate() { return interpolate_binary_search_cpp(E_target, *this); }
+```
+
+This two-step process ensures that:
+
+- The most efficient method is selected based on the data characteristics
+- Users can simply call `material.interpolate(energy_density)` without worrying about which technique to use
+
+This method:
+- Uses Double Lookup if the data is suitable (uniform or near-uniform)
+- Falls back to Binary Search for non-uniform data
+- Provides the best performance for your specific data structure without manual selection
+- Performance is optimized automatically without manual intervention
+
+## Performance Considerations
+
+Based on performance testing with 64³ cells:
+
+- **Binary Search**: Reliable but slower for large datasets
+- **Double Lookup**: Typically 2-4x faster than Binary Search for uniform data
+- **Edge Cases**: Both methods handle values outside the data range by clamping to the nearest endpoint
+
+The test results demonstrate that Double Lookup consistently outperforms Binary Search when the data is suitable, 
+making it the preferred method for uniform temperature distributions.
+
+## Implementation Details
+
+The interpolation methods are implemented in C++ for maximum performance:
+
+- The `InterpolationArrayContainer` class generates optimized C++ code with the interpolation arrays and methods
+- The generated code includes both interpolation methods and the automatic selector
+- The implementation handles edge cases gracefully, returning boundary values for out-of-range queries
+
+## Edge Case Handling
+
+Both interpolation methods include robust handling of edge cases:
+
+- **Out-of-Range Values**: Both methods clamp values outside the defined range to the nearest endpoint
+
+- **Extreme Values**: As shown in the stress tests, both methods handle extremely large and small values correctly
+
+- **Monotonicity**: The code validates that both arrays are monotonic and that energy increases with temperature
+
+The test file demonstrates this with specific test cases:
+
+```c++
+// Test with extremely large values
+constexpr double large_E = 1.0e20;
+const double result_large_bs = bs_tests.interpolateBS(large_E);
+const double result_large_dl = dl_tests.interpolateDL(large_E);
+
+// Test with extremely small values
+constexpr double small_E = 1.0e-20;
+const double result_small_bs = bs_tests.interpolateBS(small_E);
+const double result_small_dl = dl_tests.interpolateDL(small_E);
+```
+
+## When to Use Each Method
+
+- **Binary Search**: Use for general-purpose interpolation or when data points are irregularly spaced
+- **Double Lookup**: Use when data points are uniformly or nearly uniformly spaced
+- **Automatic Selection**: Use in most cases to get the best performance without manual selection
+
+For most applications, the automatic `interpolate()` method provides the best balance of performance and reliability.

docs/explanation/material_properties.md

+# Material Properties in pymatlib
+
+This document explains the conceptual framework behind temperature-dependent material properties in pymatlib,
+how they are represented internally, and the mathematical models used for computed properties.
+
+## Conceptual Framework
+
+Material properties in pymatlib are designed around these key principles:
+
+1. **Temperature Dependence**: Most material properties vary with temperature, especially during phase transitions
+2. **Symbolic Representation**: Properties are represented as symbolic expressions for mathematical manipulation
+3. **Flexible Definition**: Properties can be defined through various methods (constants, data points, files, or computation)
+4. **Physical Consistency**: Computed properties follow established physical relationships
+
+## Internal Representation
+
+### MaterialProperty Class
+
+At the core of pymatlib's property system is the `MaterialProperty` class, which contains:
+
+- A symbolic expression (`expr`) representing the property
+- A list of assignments needed to evaluate the expression
+- Methods to evaluate the property at specific temperatures
+
+```python
+@dataclass
+class MaterialProperty:
+    expr: sp.Expr
+    assignments: List[Assignment] = field(default_factory=list)
+    
+    def evalf(self, symbol: sp.Symbol, temperature: Union[float, ArrayTypes]) -> Union[float, np.ndarray]:
+        # Evaluates the property at the given temperature
+```
+
+### Assignment Class
+
+The `Assignment` class represents intermediate calculations needed for property evaluations:
+
+```python
+@dataclass
+class Assignment:
+    lhs: sp.Symbol
+    rhs: Union[tuple, sp.Expr]
+    lhs_type: str
+```
+
+## Property Definition Methods
+
+pymatlib supports multiple ways to define material properties:
+
+1. Constant Values
+
+Properties that don't vary with temperature are defined as simple numeric values:
+
+```yaml
+thermal_expansion_coefficient: 16.3e-6
+```
+
+Internally, these are converted to constant symbolic expressions.
+
+2. Interpolated Values
+
+For properties that vary with temperature, pymatlib supports interpolation between data points:
+
+```yaml
+heat_conductivity:
+    key: [1200, 1800, 2200, 2400]  # Temperatures in Kelvin
+    val: [25, 30, 33, 35]          # Property values
+```
+
+Internally, these are represented as piecewise functions that perform linear interpolation between points.
+
+3. File-Based Properties
+   
+Properties can be loaded from external data files:
+
+```yaml
+density:
+    file: ./material_data.xlsx
+    temp_col: T (K)
+    prop_col: Density (kg/(m)^3)
+```
+
+The data is loaded and converted to an interpolated function similar to key-value pairs.
+
+4. Computed Properties
+
+Some properties can be derived from others using physical relationships:
+
+```yaml
+thermal_diffusivity: compute  # k/(ρ*cp)
+```
+
+## Computed Models
+
+pymatlib implements several physical models for computing properties:
+
+### Density by Thermal Expansion
+
+```text
+ρ(T) = ρ₀ / (1 + tec * (T - T₀))³
+```
+
+Where:
+- ρ₀ is the base density
+- T₀ is the base temperature
+- tec is the thermal expansion coefficient
+
+### Thermal Diffusivity
+
+```text
+α(T) = k(T) / (ρ(T) * cp(T))
+```
+
+Where:
+- k(T) is the thermal conductivity
+- ρ(T) is the density
+- cp(T) is the specific heat capacity
+
+### Energy Density
+
+pymatlib supports multiple models for energy density:
+
+1. **Standard Model**
+
+```text
+E(T) = ρ(T) * (cp(T) * T + L)
+```
+
+2. **Enthalpy-Based Model**
+
+```text
+E(T) = ρ(T) * (h(T) + L)
+```
+
+3. **Total Enthalpy Model**
+```text
+E(T) = ρ(T) * h(T)
+```
+
+Where:
+- ρ(T) is the density
+- cp(T) is the specific heat capacity
+- h(T) is the specific enthalpy
+- L is the latent heat of fusion
+
+## Interpolation Between Data Points
+
+For properties defined through key-value pairs or files, pymatlib performs linear interpolation between data points:
+
+1. For a temperature T between two known points T₁ and T₂:
+
+```text
+property(T) = property(T₁) + (property(T₂) - property(T₁)) * (T - T₁) / (T₂ - T₁)
+```
+
+2. For temperatures outside the defined range, the property value is clamped to the nearest endpoint.
+
+## Temperature Arrays for EnergyDensity
+
+When using computed energy density, you must specify a temperature array:
+
+```text
+energy_density_temperature_array: (300, 3000, 541)  # 541 points from 300K to 3000K
+```
+
+This array is used to pre-compute energy density values for efficient interpolation during simulations.
+
+### Best Practices
+
+1. **Use Consistent Units**: All properties should use SI units (m, s, kg, K, etc.)
+2. **Cover the Full Temperature Range**: Ensure properties are defined across the entire temperature range of your simulation
+3. **Add Extra Points Around Phase Transitions**: For accurate modeling, use more data points around phase transitions
+4. **Validate Against Experimental Data**: When possible, compare property values with experimental measurements
+5. **Document Property Sources**: Keep track of where property data comes from for reproducibility

docs/how-to/define_materials.md

+# Defining Custom Material Properties
+
+This guide explains how to define custom material properties in pymatlib using different methods.
+
+## YAML Configuration Options
+
+pymatlib supports several ways to define material properties in YAML files, following SI units (m, s, kg, A, V, K, etc.).
+
+### 1. Constant Value
+
+For properties that don't vary with temperature:
+
+```yaml
+properties:
+    thermal_expansion_coefficient: 16.3e-6
+```
+
+### 2. Key-Value Pairs for Interpolation
+
+For properties that vary with temperature:
+
+```yaml
+properties:
+    # Using explicit temperature list
+    heat_conductivity:
+        key: [1200, 1800, 2200, 2400]  # Temperatures in Kelvin
+        val: [25, 30, 33, 35]  # Property values
+    
+    # Using references to defined temperatures
+    latent_heat_of_fusion:
+        key: [solidus_temperature, liquidus_temperature]
+        val: [171401, 0]
+    
+    # Using tuple for temperature generation
+    heat_capacity:
+        key: (1000, 200)  # Start at 1000K and increment by 200K for each value in val
+        # Generates: [1000, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2600, 2800, 3000, 3200]
+        val: [580, 590, 600, 600, 600, 610, 620, 630, 660, 700, 750, 750]
+    
+    # Using tuple with negative increment
+    density:
+        key: (1735.00, -5)  # Start at 1735.00K and decrement by 5K for each value in val
+        # Generates: [1735.00, 1730.00, 1725.00, 1720.00, 1715.00, 1710.00, 1705.00, 1700.00, 1695.00, 1690.00]
+        val: [7037.470, 7060.150, 7088.800, 7110.460, 7127.680, 7141.620, 7156.800, 7172.590, 7184.010, 7192.780]
+```
+
+### 3. Loading from External Files
+
+For properties defined in spreadsheets:
+
+```yaml
+properties:
+    # Simple format (first column = temperature, second column = property)
+    heat_capacity: ./heat_capacity_data.txt
+    
+    # Advanced format (specify columns)
+    density:
+      file: ./304L_data.xlsx
+      temp_col: T (K)
+      prop_col: Density (kg/(m)^3)
+```
+
+Supported file formats include .txt (space/tab separated), .csv, and .xlsx.
+
+### 4. Energy density temperature arrays
+
+For properties that need to be evaluated at specific temperature points:
+
+```yaml
+properties:
+    # Using count (integer)
+    energy_density_temperature_array: (300, 3000, 541)  # 541 evenly spaced points
+    # OR
+    # Using step size (float)
+    energy_density_temperature_array: (300, 3000, 5.0)  # From 300K to 3000K in steps of 5K
+    # OR
+    # Descending order
+    energy_density_temperature_array: (3000, 300, -5.0)  # From 3000K to 300K in steps of -5K
+```
+
+### 5. Computed Properties
+
+For properties that can be derived from others:
+
+```yaml
+properties:
+    # Simple format for density
+    density: compute
+
+    # Simple format for thermal_diffusivity
+    thermal_diffusivity: compute # Will be calculated from k/(ρ*cp)
+
+    # Simple format for energy_density
+    energy_density: compute  # Uses default model: ρ(T) * (cp(T) * T + L)
+    # OR
+    # Advanced format with explicit model selection for energy_density
+    energy_density:
+        compute: enthalpy_based  # Uses model: ρ(T) * (h(T) + L)
+    # OR
+    energy_density:
+        compute: total_enthalpy  # Uses model: ρ(T) * h(T)
+```
+
+### The equations used for computed properties are:
+- Density by thermal expansion: ρ(T) = ρ₀ / (1 + tec * (T - T₀))³
+- - Required: base_temperature, base_density, thermal_expansion_coefficient
+
+- Thermal diffusivity: α(T) = k(T) / (ρ(T) * cp(T))
+- - Required: heat_conductivity, density, heat_capacity
+
+- Energy density (standard model): ρ(T) * (cp(T) * T + L)
+- - Required: density, heat_capacity, latent_heat_of_fusion
+
+- Energy density (enthalpy_based): ρ(T) * (h(T) + L)
+- - Required: density, specific_enthalpy, latent_heat_of_fusion
+
+- Energy density (total_enthalpy): ρ(T) * h(T)
+- - Required: density, specific_enthalpy
+
+## Creating a Complete Alloy Definition
+
+Here's a complete example for stainless steel [SS304L](https://i10git.cs.fau.de/rahil.doshi/pymatlib/-/blob/master/src/pymatlib/data/alloys/SS304L/SS304L.yaml?ref_type=heads):
+
+```yaml
+name: SS304L
+
+composition:
+    C: 0.0002
+    Si: 0.0041
+    Mn: 0.016
+    P: 0.00028
+    S: 0.00002
+    Cr: 0.1909
+    N: 0.00095
+    Ni: 0.0806
+    Fe: 0.70695
+
+solidus_temperature: 1605
+liquidus_temperature: 1735
+
+properties:
+    energy_density:
+      compute: total_enthalpy
+    energy_density_temperature_array: (300, 3000, 541)
+
+    base_temperature: 2273.
+    base_density: 6.591878918e3
+
+    density:
+      file: ./304L_data.xlsx
+      temp_col: T (K)
+      prop_col: Density (kg/(m)^3)
+
+    heat_conductivity:
+      file: ./304L_data.xlsx
+      temp_col: T (K)
+      prop_col: Thermal conductivity (W/(m*K))
+
+      heat_capacity:
+        file: ./304L_data.xlsx
+        temp_col: T (K)
+        prop_col: Specific heat (J/(Kg K))
+
+      thermal_expansion_coefficient: 16.3e-6
+    
+      specific_enthalpy:
+        file: ./304L_data.xlsx
+        temp_col: T (K)
+        prop_col: Enthalpy (J/kg)
+    
+      latent_heat_of_fusion:
+        key: [solidus_temperature, liquidus_temperature]
+        val: [171401, 0]
+    
+      thermal_diffusivity: compute
+```
+
+## Best Practices
+
+- Use consistent units throughout your definitions
+- Document the expected units for each property
+- For temperature-dependent properties, cover the full range of temperatures you expect in your simulation
+- Validate your property data against experimental values when possible
+- Use computed properties only when the relationship is well-established
+- All numerical values must use period (.) as decimal separator, not comma
+- Interpolation between data points is performed automatically for file-based and key-val properties
+
+## Important Notes
+
+- If a specific property is defined in multiple ways or multiple times, the parser will throw an error
+- If required dependencies for computed properties are missing, an error will be raised
+- Properties will be computed in the correct order regardless of their position in the file
+- To retrieve temperature from energy_density, use the default "interpolate" method from within the generated class from InterpolationArrayContainer named after the alloy

docs/how-to/energy_temperature_conversion.md

+# Converting Between Energy Density and Temperature
+
+This guide explains how to perform bilateral conversions between energy density and temperature in pymatlib.
+
+## Why Energy-Temperature Conversion Matters
+
+In many material simulations, particularly those involving heat transfer and phase changes, you need to convert between:
+
+- **Temperature**: The conventional measure of thermal state
+- **Energy Density**: The amount of thermal energy per unit volume
+
+This conversion is essential because:
+- Energy is conserved in physical processes
+- Phase transitions involve latent heat where temperature remains constant
+- Many numerical methods work better with energy as the primary variable
+
+## Basic Conversion Process
+
+### 1. Generate the Interpolation Container
+
+First, create an interpolation container that stores the energy-temperature relationship:
+
+```python
+import pystencils as ps
+from pystencilssfg import SourceFileGenerator
+from pymatlib.core.yaml_parser import create_alloy_from_yaml
+from pymatlib.core.interpolators import InterpolationArrayContainer
+
+with SourceFileGenerator() as sfg:
+u = ps.fields("u: float64[2D]", layout='fzyx')
+
+# Create an alloy
+alloy = create_alloy_from_yaml("path/to/alloy.yaml", u.center())
+
+# Generate the interpolation container
+arr_container = InterpolationArrayContainer.from_material("SS304L", alloy)
+sfg.generate(arr_container)
+```
+
+This generates C++ code with the `InterpolationArrayContainer` class that contains:
+- Temperature array
+- Energy density array
+- Methods for interpolation
+
+### 2. Using the Generated Code
+
+In your C++ application, you can use the generated code:
+
+```c++
+// Energy to temperature conversion using binary search
+double energy_density = 16950000000.0;  // J/m³
+SS304L material;
+double temp = material.interpolateBS(energy_density);
+
+// Energy to temperature conversion using double lookup
+double temp2 = material.interpolateDL(energy_density);
+
+// Using the automatic method selection
+double temp3 = material.interpolate(energy_density);
+```
+
+Note that the temperature to energy density conversion is handled in Python using the evalf method:
+
+```python
+import sympy as sp
+
+# Create symbolic temperature variable
+T = sp.Symbol('T')
+
+# Temperature to energy conversion in Python
+temperature = 1500.0  # Kelvin
+energy = alloy.energy_density.evalf(T, temperature)
+
+```
+
+## Interpolation Methods
+
+pymatlib provides two interpolation methods:
+
+### Binary Search Interpolation
+
+```c++
+double temp = material.interpolateBS(energy_density);
+```
+
+- Best for non-uniform temperature distributions
+- O(log n) lookup complexity
+- Robust for any monotonically increasing data
+
+### Double Lookup Interpolation
+
+```c++
+double temp = material.interpolateDL(energy_density);
+```
+
+- Optimized for uniform temperature distributions
+- O(1) lookup complexity
+- Faster but requires pre-processing
+
+### Automatic Method Selection
+
+```c++
+double temp = material.interpolate(energy_density);
+```
+
+## Custom Interpolation Arrays
+
+You can create custom interpolation containers for specific temperature ranges:
+
+```python
+import numpy as np
+from pymatlib.core.interpolators import InterpolationArrayContainer
+
+# Create custom temperature and energy arrays
+T_array = np.linspace(300, 3000, 5)
+E_array = np.array([...]) # Your energy values
+
+# Create a custom container
+custom_container = InterpolationArrayContainer("CustomMaterial", T_array, E_array)
+sfg.generate(custom_container)
+```
+
+## Performance Considerations
+
+- **Binary Search**: Use for general-purpose interpolation
+- **Double Lookup**: Use for uniform temperature distributions
+- **Array Size**: Balance between accuracy and memory usage
+- **Density**: Add more points around phase transitions for accuracy
+
+## Complete Example
+
+Here's a complete example showing both methods:
+
+```python
+import numpy as np
+from pystencilssfg import SourceFileGenerator
+from pymatlib.core.interpolators import InterpolationArrayContainer
+
+with SourceFileGenerator() as sfg:
+    # Create temperature and energy arrays
+    T_bs = np.array([3243.15, 3248.15, 3258.15, 3278.15], dtype=np.float64)
+    E_bs = np.array([1.68e10, 1.69e10, 1.70e10, 1.71e10], dtype=np.float64)
+
+    # Binary search container
+    binary_search_container = InterpolationArrayContainer("BinarySearchTests", T_bs, E_bs)
+    sfg.generate(binary_search_container)
+    
+    # Double lookup container
+    T_eq = np.array([3243.15, 3253.15, 3263.15, 3273.15], dtype=np.float64)
+    E_neq = np.array([1.68e10, 1.69e10, 1.70e10, 1.71e10], dtype=np.float64)
+    double_lookup_container = InterpolationArrayContainer("DoubleLookupTests", T_eq, E_neq)
+    sfg.generate(double_lookup_container)
+```

docs/reference/api/alloy.md

+## API Reference
+
+### Core Classes
+
+#### Alloy
+
+A dataclass representing a material alloy.
+
+**Properties:**
+- `elements`: List of chemical elements in the alloy
+- `composition`: Dictionary mapping element symbols to their weight fractions
+- `temperature_solidus`: Solidus temperature (K)
+- `temperature_liquidus`: Liquidus temperature (K)
+- `properties`: Dictionary of material properties

docs/reference/yaml_schema.md

+# YAML Schema for Material Definition
+
+This document defines the schema for material definition YAML files in pymatlib.
+
+## Schema Overview
+
+A valid material definition must include:
+- `name`: String identifier for the material
+- `composition`: Map of element symbols to their mass fractions
+- `solidus_temperature`: Numeric value in Kelvin
+- `liquidus_temperature`: Numeric value in Kelvin
+- `properties`: Map of property names to their definitions
+
+## Top-Level Structure
+
+```yaml
+name: <string> # Required: Material name
+composition: # Required: Chemical composition
+<element>: <mass_fraction> # At least one element required
+solidus_temperature: <float> # Required: Solidus temperature in K
+liquidus_temperature: <float> # Required: Liquidus temperature in K
+properties: # Required: Material properties
+<property_name>: <definition> # At least one property required
+```
+
+## Property Definition Types
+
+Properties can be defined in four different ways:
+
+### 1. Constant Value
+
+For properties that don't vary with temperature:
+
+```yaml
+properties:
+    thermal_expansion_coefficient: 16.3e-6 # Single numeric value
+```
+
+### 2. File-Based Properties
+
+#### 2.1 Simple Format
+
+```yaml
+properties:
+    density: ./density_temperature.txt # Path to file with two columns
+```
+
+The simple format assigns the first column to temperature data (K) and the second column to the corresponding property values.
+
+#### 2.2 Advanced Format
+
+```yaml
+properties:
+    density:
+        file: ./path/to/file.xlsx
+        temp_col: Temperature # Name of the temperature column
+        prop_col: Property # Name of the property column
+```
+
+Advanced format is required when you have a file with multiple columns.
+Supported file formats: .txt (space/tab separated), .csv, and .xlsx.
+
+### 3. Key-Value Pairs
+
+For properties defined at specific temperature points:
+
+```yaml
+properties:
+    density:
+        key: [1605, 1735] # Temperature values in K
+        val: [7262.34, 7037.47] # Corresponding property values
+```
+
+Three ways to specify temperature points:
+
+#### 3.1 Explicit List
+
+```yaml
+key: [1735.00, 1730.00, 1720.00, 1715.00] # Explicit temperature list
+val: [7037.470, 7060.150, 7110.460, 7127.680] # Corresponding values
+```
+
+#### 3.2 Reference to Defined Temperatures
+
+```yaml
+key: [solidus_temperature, liquidus_temperature] # References
+val: [7262.34, 7037.47] # Corresponding values
+```
+
+#### 3.3 Tuple with Start and Increment
+
+```yaml
+key: (1735.00, -5) # Start at 1735.00K and decrement by 5K
+val: [7037.470, 7060.150, 7088.800, 7110.460, 7127.680]
+```
+
+When using a tuple for key, the generated temperature points will be:
+`[start, start+increment, start+2*increment, ...]` until matching the length of val.
+
+### 4. Computed Properties
+
+#### 4.1 Simple Format
+
+```yaml
+properties:
+    density: compute
+    thermal_diffusivity: compute
+    energy_density: compute
+```
+
+Simple format uses the default model to compute the property.
+
+#### 4.2 Advanced Format
+
+```yaml
+properties:
+    energy_density:
+        compute: enthalpy_based # Specific computation model
+```
+
+#### 4.3 Energy Density Temperature Array
+
+When energy_density is computed, you must specify the temperature array:
+
+```yaml
+properties:
+    energy_density: compute
+    energy_density_temperature_array: (300, 3000, 541) # 541 points
+```
+
+The third parameter can be:
+- An integer: Total number of points to generate
+- A float: Temperature increment/decrement between points
+
+## Computation Models and Required Properties
+
+### Density (by thermal expansion)
+- **Equation**: ρ(T) = ρ₀ / (1 + tec * (T - T₀))³
+- **Required properties**: base_temperature, base_density, thermal_expansion_coefficient
+
+### Thermal Diffusivity
+- **Equation**: α(T) = k(T) / (ρ(T) * cp(T))
+- **Required properties**: heat_conductivity, density, heat_capacity
+
+### Energy Density (standard model)
+- **Equation**: ρ(T) * (cp(T) * T + L)
+- **Required properties**: density, heat_capacity, latent_heat_of_fusion
+
+### Energy Density (enthalpy_based model)
+- **Equation**: ρ(T) * (h(T) + L)
+- **Required properties**: density, specific_enthalpy, latent_heat_of_fusion
+
+### Energy Density (total_enthalpy model)
+- **Equation**: ρ(T) * h(T)
+- **Required properties**: density, specific_enthalpy
+
+## Validation Rules
+
+1. All required top-level fields must be present
+2. Composition fractions must sum to approximately 1.0
+3. Liquidus temperature must be greater than or equal to solidus temperature
+4. Properties cannot be defined in multiple ways or multiple times
+5. Required dependencies for computed properties must be present
+6. Temperature arrays must be monotonic
+7. Energy density arrays must be monotonic with respect to temperature
+8. File paths must be valid and files must exist
+9. For key-value pairs, key and val must have the same length
+10. When using tuple notation for temperature arrays, the increment must be non-zero
+
+## Important Notes
+
+1. All numerical values must use period (.) as decimal separator, not comma
+2. Interpolation between data points is performed automatically for file-based and key-val properties
+3. Properties will be computed in the correct order regardless of their position in the file
+4. To retrieve temperature from energy_density, use the default "interpolate" method from within the generated class

docs/tutorials/first_simulation.md

+# Creating Your First Material Simulation
+
+This tutorial will guide you through creating a basic heat equation simulation using [pymatlib](https://i10git.cs.fau.de/rahil.doshi/pymatlib) and [pystencils](https://pycodegen.pages.i10git.cs.fau.de/pystencils/).
+It builds upon the existing [waLBerla tutorial for code generation](https://walberla.net/doxygen/tutorial_codegen01.html), adding material property handling with pymatlib.
+
+## Prerequisites
+
+Before starting, ensure you have:
+- Completed the [Getting Started](getting_started.md) tutorial
+- Installed pystencils and [pystencilssfg](https://pycodegen.pages.i10git.cs.fau.de/pystencils-sfg/)
+- Basic understanding of heat transfer and the heat equation
+- Created the `simple_steel.yaml` file from the [Getting Started](getting_started.md) tutorial
+
+## Overview
+
+We'll create a 2D heat equation simulation with temperature-dependent material properties.
+The main enhancement compared to the original waLBerla tutorial is that we'll use pymatlib to handle material properties that vary with temperature.
+
+The steps are:
+
+1. Define a temperature field
+2. Create an alloy with temperature-dependent properties
+3. Set up the heat equation with material properties
+4. Generate code for the simulation
+5. Run the simulation
+
+## Step 1: Define the Temperature Field
+
+First, we'll create a temperature field using pystencils:
+
+```python
+import pystencils as ps
+from pystencilssfg import SourceFileGenerator
+
+with SourceFileGenerator() as sfg:
+    data_type = "float64"
+
+    # Define temperature fields and output field for thermal diffusivity
+    u, u_tmp = ps.fields(f"u, u_tmp: {data_type}[2D]", layout='fzyx')
+    thermal_diffusivity_out = ps.fields(f"thermal_diffusivity_out: {data_type}[2D]", layout='fzyx')
+```
+
+## Step 2: Set Up the Heat Equation
+
+Now we'll set up the heat equation with temperature-dependent material properties:
+
+```python
+import sympy as sp
+
+# Create symbolic variables for the equation
+dx = sp.Symbol("dx")
+dt = sp.Symbol("dt")
+thermal_diffusivity = sp.Symbol("thermal_diffusivity")
+
+# Define the heat equation using finite differences
+heat_pde = ps.fd.transient(u) - thermal_diffusivity * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
+
+# Discretize the PDE
+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()
+```
+
+## Step 3: Create an Alloy
+
+This is where pymatlib enhances the original tutorial.
+We'll load the alloy from a YAML file:
+
+```python
+from pymatlib.core.yaml_parser import create_alloy_from_yaml
+
+# Load the SimpleSteel alloy definition we created earlier
+simple_steel_alloy = create_alloy_from_yaml("path/to/simple_steel.yaml", u.center())
+```
+
+The second parameter `u.center()` links the alloy to the temperature field, making material properties dependent on the local temperature.
+
+## Step 4: Convert Material Property Assignments
+
+We need to convert the material property assignments to a format pystencils can understand:
+
+```python
+from pymatlib.core.assignment_converter import assignment_converter
+
+# Access the thermal diffusivity property directly from the material
+# The property is already defined in the YAML file or computed automatically
+
+# Convert material property assignments to pystencils format
+subexp, subs = assignment_converter(simple_steel_alloy.thermal_diffusivity.assignments)
+
+# Add the thermal diffusivity expression to subexpressions
+subexp.append(ps.Assignment(thermal_diffusivity, simple_steel_alloy.thermal_diffusivity.expr))
+```
+
+The `assignment_converter` function transforms pymatlib's internal representation of material properties into pystencils assignments. 
+It handles type conversions and creates properly typed symbols for the code generation.
+
+## Step 5: Generate Code for the Simulation
+
+Finally, we'll create the assignment collection and generate optimized C++ code for the simulation:
+
+```python
+from sfg_walberla import Sweep
+from pymatlib.core.interpolators import InterpolationArrayContainer
+
+# Create assignment collection with subexpressions and main assignments
+ac = ps.AssignmentCollection(
+    subexpressions=subexp,  # Include the subexpressions from pymatlib
+    main_assignments=[
+        ps.Assignment(u_tmp.center(), heat_pde_discretized),
+        ps.Assignment(thermal_diffusivity_out.center(), thermal_diffusivity)
+    ])
+
+# Generate the sweep
+sweep = Sweep("HeatEquationKernelWithMaterial", ac)
+sfg.generate(sweep)
+
+# Create the container for energy-temperature conversion
+# The name "SimpleSteel" is just an identifier for the generated class
+arr_container = InterpolationArrayContainer.from_material("SimpleSteel", simple_steel_alloy)
+sfg.generate(arr_container)
+```
+
+## Step 6: Create the Interpolation Container
+
+For energy-temperature conversion, create an interpolation container:
+
+```python
+from pymatlib.core.interpolators import InterpolationArrayContainer
+
+# Create the container for energy-temperature conversion
+arr_container = InterpolationArrayContainer.from_material("SS304L", alloy)
+sfg.generate(arr_container)
+```
+
+## Complete Example
+
+Here's the complete example:
+
+```python
+import sympy as sp
+import pystencils as ps
+from pystencilssfg import SourceFileGenerator
+from sfg_walberla import Sweep
+from pymatlib.core.assignment_converter import assignment_converter
+from pymatlib.core.interpolators import InterpolationArrayContainer
+from pymatlib.core.yaml_parser import create_alloy_from_yaml
+
+with SourceFileGenerator() as sfg:
+    data_type = "float64"
+
+    # Define fields
+    u, u_tmp = ps.fields(f"u, u_tmp: {data_type}[2D]", layout='fzyx')
+    thermal_diffusivity = sp.Symbol("thermal_diffusivity")
+    thermal_diffusivity_out = ps.fields(f"thermal_diffusivity_out: {data_type}[2D]", layout='fzyx')
+    dx, dt = sp.Symbol("dx"), sp.Symbol("dt")
+
+    # Define heat equation
+    heat_pde = ps.fd.transient(u) - thermal_diffusivity * (ps.fd.diff(u, 0, 0) + ps.fd.diff(u, 1, 1))
+
+    # Discretize the PDE
+    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()
+
+    # Create alloy with temperature-dependent properties
+    mat = create_alloy_from_yaml("simple_steel.yaml", u.center())
+
+    # Convert material property assignments to pystencils format
+    subexp, subs = assignment_converter(mat.thermal_diffusivity.assignments)
+    subexp.append(ps.Assignment(thermal_diffusivity, mat.thermal_diffusivity.expr))
+
+    # Create assignment collection with the converted subexpressions
+    ac = ps.AssignmentCollection(
+        subexpressions=subexp,
+        main_assignments=[
+            ps.Assignment(u_tmp.center(), heat_pde_discretized),
+            ps.Assignment(thermal_diffusivity_out.center(), thermal_diffusivity)
+        ])
+
+    # Generate the sweep
+    sweep = Sweep("HeatEquationKernelWithMaterial", ac)
+    sfg.generate(sweep)
+
+    # Generate interpolation container for energy-temperature conversion
+    arr_container = InterpolationArrayContainer.from_material("SimpleSteel", mat)
+    sfg.generate(arr_container)
+```
+
+## Understanding the Assignment Converter
+
+The `assignment_converter` function is a key component of pymatlib that bridges the gap between symbolic material property representations and pystencils code generation. It:
+
+1. Takes a list of `Assignment` objects from pymatlib
+2. Converts them to pystencils-compatible assignments with proper typing
+3. Returns both the converted assignments and a mapping of symbols
+
+This allows material properties defined in YAML files to be seamlessly integrated into pystencils code generation, enabling temperature-dependent properties in simulations.
+
+The main ways pymatlib enhances the simulation are:
+
+- Material property handling: Properties are defined in YAML and accessed via the material object
+- Assignment conversion: The assignment_converter transforms material property assignments to pystencils format
+- Energy-temperature conversion: The InterpolationArrayContainer enables efficient conversion between energy density and temperature
+
+## Next Steps
+
+Now that you've created your first simulation with temperature-dependent material properties, you can:
+- Learn about [energy-temperature conversion](../how-to/energy_temperature_conversion.md)
+- Explore more complex [material properties](../explanation/material_properties.md)
+- Understand the [API reference](../reference/api/alloy.md) for advanced usage
+- Check the original [waLBerla tutorial](https://walberla.net/doxygen/tutorial_codegen01.html) for details on running the simulation

docs/tutorials/getting_started.md

+# Getting Started with pymatlib
+
+This tutorial will guide you through the basics of using pymatlib to model material properties for simulations.
+
+## Prerequisites
+
+Before starting, ensure you have:
+
+- Python 3.10 or newer
+- Basic knowledge of material properties
+- Basic familiarity with Python
+
+## Installation
+
+Install pymatlib using pip:
+
+```bash
+pip install "git+https://i10git.cs.fau.de/rahil.doshi/pymatlib.git"
+```
+
+
+For development, clone the repository and install in development mode:
+
+```bash
+git clone https://i10git.cs.fau.de/rahil.doshi/pymatlib.git
+cd pymatlib
+pip install -e .
+```
+
+## Basic Concepts
+
+pymatlib organizes material data around a few key concepts:
+
+1. **Alloys**: Compositions of multiple elements with specific properties
+2. **Material Properties**: Physical characteristics that can vary with temperature
+3. **Interpolation**: Methods to estimate property values between known data points
+
+## Your First Alloy
+
+Let's create a simple alloy definition:
+
+1. Create a file named `simple_steel.yaml` with the following content:
+```python
+name: SimpleSteel
+
+composition:
+    Fe: 0.98
+    C: 0.02
+
+solidus_temperature: 1450
+liquidus_temperature: 1520
+
+properties:
+    density:
+        key:[300, 800, 1300, 1800]
+        val:[7850, 7800, 7750, 7700]
+
+    heat_conductivity:
+        key: [300, 800, 1300, 1800]
+        val: [18.5, 25, 32, 36.5]
+
+    heat_capacity:
+        key: [300, 800, 1300, 1800]
+        val: [450, 500, 550, 600]
+
+    thermal_diffusivity: compute    
+```
+
+2. Load the alloy in Python:
+```python
+from pymatlib.core.yaml_parser import create_alloy_from_yaml
+
+# Load the alloy definition
+alloy = create_alloy_from_yaml("simple_steel.yaml")
+
+# Print basic information
+print(f"Alloy: {alloy.name}")
+print(f"Composition: {alloy.composition}")
+print(f"Temperature range: {alloy.temperature_solidus}K - {alloy.temperature_liquidus}K")
+
+# Get property values at specific temperatures
+temp = 500 # Kelvin
+density = alloy.get_property("density", temp)
+conductivity = alloy.get_property("heat_conductivity", temp)
+capacity = alloy.get_property("heat_capacity", temp)
+
+print(f"At {temp}K:")
+print(f" Density: {density} kg/m³")
+print(f" Thermal Conductivity: {conductivity} W/(m·K)")
+print(f" Heat Capacity: {capacity} J/(kg·K)")
+```
+
+## Next Steps
+
+Now that you've created your first alloy, you can:
+
+- Learn how to [create your first simulation](first_simulation.md)
+- Explore [defining custom material properties](../how-to/define_materials.md)
+- Understand [interpolation methods](../explanation/interpolation_methods.md)