ENH: Workspace API

[richardotis]

• Adds a new Workspace object and associated imperative API. This object now underlies the equilibrium function, which remains for backwards compatibility. A Workspace object can be mutated after creation. The Workspace.get() function accepts ComputableProperty objects as input and returns arrays; the Workspace.plot() function works similarly, but returns a matplotlib figure. Fixes Develop a wrapper API for plotting properties from calculate and equilibrium #154
• Adds a new PhaseRecordFactory object and associated connections to PhaseRecord, which enables deferred compilation of symbolically-defined model properties until they are called in the PhaseRecord API. The PhaseRecord object also gains prop and prop_grad functions for computing properties and property gradients of arbitrary symbolic attributes of Model instances.
• Creates a "Computable Property Framework" (CPF) for calculating properties of the system, individual phases, or some combination of these. It also provides a framework for defining custom properties.
• Adds support for "dot" derivatives (constrained total derivatives) of thermodynamic properties (fixes Heat capacities calculated by Model can sometimes be incorrect #261 )
• Adds support for T0 (tzero) calculation
• Adds support for driving force calculation, including nucleation driving forces, via the IsolatedPhase and DormantPhase meta-properties. Meta-properties are objects that know how to create new properties (as defined by the CPF). In this case, these meta-properties know how to start sub-calculations . The current way of finding starting points for these calculations is not ideal and prone to failure (basically reuses CompositionSet objects from the workspace equilibrium calculation), but it works well in the typical case.
• Adds physical units support via pint. All ComputableProperty objects have "implementation units" and "display units." Implementation units are what are assumed by the underlying numerical computation (which is done without units, for performance). The display units are what the user gets when Workspace.get or Workspace.plot are called. ComputableProperty.__getitem__ allows the display units to be changed for individual calls to Workspace.get, Workspace.plot, or when setting Workspace.conditions.
• Also regarding physical units, a special unit conversion context is defined to allow conversions between per-mole properties and per-mass properties. There's room to add per-volume here as well, but that may not make it into this PR.
• Adds support for specifying linear combinations and ratios of mole fractions as a condition.
• Adds support for specifying weight fractions as a condition.
• Associated tests for these new features.

Documentation

• Adds several examples for the Computable Property API, including driving force, T0, and isolated energy surfaces. A demo of export-to-DataFrame is also shown using the new PandasRenderer.
• Other examples have been updated to use the new APIs. Examples using the legacy APIs continue to be supported, but have been moved to a new "Advanced Examples" section. These sections will continue to be expanded and reorganized to accommodate our evolving best practices.

• If any dependencies have changed, the changes are reflected in the
  • setup.py (runtime requirements)

[bocklund]

I'm just starting to use this a bit, I'm wondering if I missed if there's a hook to give units for ModelComputedProperty that are not derived from the energy. Specifically, I've been looking at VM, V0, etc

[richardotis]

I'm just starting to use this a bit, I'm wondering if I missed if there's a hook to give units for ModelComputedProperty that are not derived from the energy. Specifically, I've been looking at VM, V0, etc

The way it works at the moment is somewhat hacky. You can assign a variable to this module: https://github.com/richardotis/pycalphad/blob/a0541fe69317348d0627ab72440b00c84a3891a0/pycalphad/property_framework/units.py#L22-L36

[bocklund]

I came across a bug for ternaries when using the cpf for mole fractions specified to be greater than one. Here's some code to reproduce:

@select_database("alnipt.tdb")
def test_jansson_derivative_with_invalid_mass_conditions(load_database):
    """
    CPF values including Jansson derivatives computed for conditions that are invalid should produce NaN.
    """
    # change to load_database() for the actual test
    dbf = load_database()
    wks = Workspace(dbf, ["AL", "NI", "PT"], ["LIQUID"], {v.T: 298.15, v.P: 101325, v.N: 1, v.X("AL"): 0.6, v.X("PT"): 0.6})
    T = wks.get("T")
    assert np.isnan(T)
    GM = wks.get("GM")
    assert np.isnan(GM)
    dGM_dT = wks.get("GM.T")
    assert np.isnan(dGM_dT)

# The following code can be used to test interactively
import numpy as np
from pycalphad import variables as v
from pycalphad.core.workspace import Workspace
from pycalphad import Database
from importlib_resources import files
import pycalphad.tests.databases
dbf = Database(str(files(pycalphad.tests.databases).joinpath("alnipt.tdb")))
wks = Workspace(dbf, ["AL", "NI", "PT"], ["LIQUID"], {v.T: 298.15, v.P: 101325, v.N: 1, v.X("AL"): 0.6, v.X("PT"): 0.6})
T = wks.get("T")  # currently NaN
print(T)
GM = wks.get("GM")  # currently 0
print(GM)
dGM_dT = wks.get("GM.T")  # raises exception
print(dGM_dT)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[23], line 30
     28 GM = wks.get("GM")  # currently 0
     29 print(GM)
---> 30 dGM_dT = wks.get("GM.T")  # raises exception
     31 print(dGM_dT)

File ~/src/pycalphad-workspace-development/pycalphad/pycalphad/core/workspace.py:425, in Workspace.get(self, values_only, *args)
    423         if results.get(arg, None) is None:
    424             results[arg] = np.zeros((arr_size,) + arg.shape)
--> 425         results[arg][local_index, ...] = Q_(arg.compute_property(composition_sets, cur_conds, chemical_potentials),
    426                                             prop_implementation_units).to(prop_display_units, context).magnitude
    427     local_index += 1
    429 for arg in args:

File ~/src/pycalphad-workspace-development/pycalphad/pycalphad/property_framework/computed_property.py:288, in DotDerivativeComputedProperty.compute_property(self, compsets, cur_conds, chemical_potentials)
    286 solver = Solver()
    287 print(compsets)
--> 288 spec = solver.get_system_spec(compsets, cur_conds)
    289 state = spec.get_new_state(compsets)
    290 state.chemical_potentials[:] = chemical_potentials

File ~/src/pycalphad-workspace-development/pycalphad/pycalphad/core/solver.py:54, in Solver.get_system_spec(self, composition_sets, conditions)
     52 from pycalphad.variables import ChemicalPotential, MoleFraction, SiteFraction
     53 compsets = composition_sets
---> 54 state_variables = compsets[0].phase_record.state_variables
     55 nonvacant_elements = compsets[0].phase_record.nonvacant_elements
     56 num_statevars = len(state_variables)

IndexError: list index out of range

It doesn't currently seem to affect binaries, i.e.

X_AL = 2.0
wks = Workspace(dbf, ["AL", "NI"], ["LIQUID"], {v.T: 298.15, v.P: 101325, v.N: 1, v.X("AL"): X_AL})

seems to work for any value of X_AL and you get results as if X_AL = 1.0.

It also breaks for tuple conditions as if they were passed as scalars, for example

wks = Workspace(dbf, ["AL", "NI", "PT"], ["LIQUID"], {v.T: 298.15, v.P: 101325, v.N: 1, v.X("AL"): (0, 1, 0.01), v.X("PT"): 0.98})

works for the first 2 conditions [0.00 0.02 0.98], [0.01 0.01 0.98], but fails on the third

[bocklund]

The test assumes that the value of "GM" should also be NaN for the invalid conditions, while it's currently zero. NaN is more consistent with the existing pre-workspace behavior, but if there's a good reason for it to be zero, that might be fine