Source code for pycalphad.property_framework.units

import pint
import numpy as np
import numpy.typing as npt
from typing import TYPE_CHECKING, Sequence
from pycalphad.core.composition_set import CompositionSet
if TYPE_CHECKING:
    from pycalphad.property_framework import ComputableProperty

ureg = pint.UnitRegistry(
    preprocessors=[lambda s: s.replace('%', ' percent ')],
    # Suppress warnings when redefining 'percent', '%', and 'ppm' units.
    # We intentionally redefine these relative to our custom 'fraction' unit.
    on_redefinition='ignore',
)
ureg.define('atom = 1/avogadro_number * mol')
ureg.define('fraction = []')
ureg.define('percent = 1e-2 fraction = %')
ureg.define('ppm = 1e-6 fraction')
pint.set_application_registry(ureg)
Q_ = ureg.Quantity
DimensionalityError = pint.DimensionalityError



def as_quantity(prop: "ComputableProperty", qt: npt.ArrayLike):
    if not isinstance(qt, Q_):
        return Q_(qt, prop.display_units)
    else:
        return qt


energy_implementation_units = GM_implementation_units = 'J / mol'
energy_display_units = GM_display_units = 'J / mol'
energy_display_name = GM_display_name = 'Gibbs Energy'
G_implementation_units = 'J'
G_display_units = 'J'
G_display_name = 'Gibbs Energy'
enthalpy_implementation_units = HM_implementation_units = GM_implementation_units
enthalpy_display_units = HM_display_units = GM_display_units
enthalpy_display_name = HM_display_name = 'Enthalpy'
H_implementation_units = 'J'
H_display_units = 'J'
H_display_name = 'Enthalpy'
entropy_implementation_units = SM_implementation_units = 'J / mol / K'
entropy_display_units = SM_display_units = 'J / mol / K'
entropy_display_name = SM_display_name = 'Entropy'
cpm_implementation_units = CPM_implementation_units = 'J / mol / K'
cpm_display_units = CPM_display_units = 'J / mol / K'
cpm_display_name = CPM_display_name = 'Heat Capacity'

def _conversions_per_formula_unit(compset):
    components = compset.phase_record.nonvacant_elements
    num_components = len(components)
    moles_per_fu = np.zeros((num_components,1))
    for comp_idx in range(num_components):
        compset.phase_record.formulamole_obj(moles_per_fu[comp_idx, :], compset.dof, comp_idx)
    # now we have 'moles per formula unit'
    # need to convert by adding molecular weight of each element
    grams_per_mol = np.array(compset.phase_record.molar_masses, dtype='float')
    grams_per_fu = np.dot(grams_per_mol, moles_per_fu)
    return moles_per_fu.sum(), grams_per_fu



def unit_conversion_context(compsets, prop):
    context = pint.Context()
    # these will be something/mol by convention
    # XXX: This is a very rough check
    if not ('/ mol' in str(prop.implementation_units)):
        return context
    implementation_units = (ureg.Unit(prop.implementation_units) * ureg.Unit('mol'))
    molar_weight = 0.0 # g/mol-atom
    for compset in compsets:
        if compset.NP > 0:
            moles_per_fu, grams_per_fu = _conversions_per_formula_unit(compset)
            grams_per_mol_atoms = (compset.NP / moles_per_fu) * grams_per_fu
            molar_weight += grams_per_mol_atoms
    molar_weight = Q_(molar_weight, 'g/mol')
    per_moles = ureg.get_dimensionality(ureg.Unit('{} / mol'.format(implementation_units)))
    per_mass = ureg.get_dimensionality(ureg.Unit('{} / g'.format(implementation_units)))

    context.add_transformation(
        per_moles,
        per_mass,
        lambda ureg, x: np.true_divide(x, molar_weight).to_reduced_units(),
    )
    context.add_transformation(
        per_mass,
        per_moles,
        lambda ureg, x: (x * molar_weight).to_reduced_units()
    )

    return context



def _composition_sets_for_unit_conversion(comp_sets: Sequence[CompositionSet], prop: "ComputableProperty") -> Sequence[CompositionSet]:
    """Select composition sets matching a phase-specific property for unit conversion context."""
    phase_name = getattr(prop, 'phase_name', None)
    if phase_name is None or phase_name == '*':
        return comp_sets

    tokens = phase_name.split('#')
    phase_name = tokens[0]
    multiplicity = int(tokens[1]) if len(tokens) > 1 else 1
    multiplicity_seen = 0
    for comp_set in comp_sets:
        if comp_set.phase_record.phase_name == phase_name:
            multiplicity_seen += 1
            if multiplicity_seen == multiplicity:
                return [comp_set]
    return comp_sets




def to_display_units(value: npt.ArrayLike, comp_sets: Sequence[CompositionSet], prop: "ComputableProperty") -> npt.ArrayLike:
    """Convert a property value from implementation units to display units.

    Filters composition sets to match phase-specific properties before
    building the molar mass context for per-mole to per-mass conversions.
    """
    implementation_units = ureg.Unit(getattr(prop, 'implementation_units', ''))
    display_units = ureg.Unit(getattr(prop, 'display_units', ''))
    context = unit_conversion_context(_composition_sets_for_unit_conversion(comp_sets, prop), prop)
    return Q_(value, implementation_units).to(display_units, context).magnitude