Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: variables.MassFraction and variables.Composition object implementations #254

Merged
merged 16 commits into from Feb 6, 2020

Conversation

@bocklund
Copy link
Collaborator

bocklund commented Dec 11, 2019

Opening this PR for feedback.

Current status

I sketched out example implementations of MassFraction and Composition objects (the old Composition object was refactored to MoleFraction).

MassFraction objects are basically sister objects to MoleFraction, just using W instead of X.

Composition objects are convenience objects that would describe a complete composition for a system as would be used in a calculation. The key ideas of a Composition are

  1. All Compositions are single point compositions.
    a. Any mixing, broadcasting, or paths between Composition objects are not Compositions themselves, but simply dicts mapping v.X or v.W to a float, tuple, or array, as in the equilibrium API.
    b. The Composition object could be a platform to create such dictionaries, but are fundamentally not those things.
  2. Can be created from mass or mole fractions and easily convert between them
  3. Dependent components can be changed
  4. Equality convenience checks between two Compositions
  5. Two compositions can be mixed (to create an new point Composition by Composition.mix) and linearly interpolated (to create a dict of {v.X/v.W: [list of values]} by Composition.interpolate_path)

The usual way I imagine a user would use this would be:

from pycalphad import Database, equilibrium, variables as v
dbf = Database('multicomponent.tdb')
comps = ['FE', 'CR', 'NI', 'MO']
phases = list(dbf.phases.keys())
ss316 = v.Composition(dbf, {v.W('CR'): 0.17, v.W('NI'): 0.12, v.W('MO'): 0.025}, 'FE')
conds = {v.P: 101325, v.T: (300, 2000, 20), v.N: 1, **ss316.mole_fractions}
eq = equilibrium(dbf, comps, phases, conds)

Some things that are still open:

  1. A v.X/v.W of a non-pure element species will currently break things. Does it make sense to/should species support be implemented?
  2. Compositions are currently instantiated by either passing a dictionary of masses or a Database (where the mass dictionary will be created from the element refdata). The main API I imagine that users will want to is Composition(dbf, {v.X('A'): ...}, 'B'). In principle Composition objects should be independent of any Database, besides a user might not want to load a Database (or might not have a database with the elements they want to manipulate, or might use fictitious elements). Should we provide a hardcoded mass dictionary for pure elements as the default and allow databases or custom mass dicts to be input as an optional argument?
  3. Should the conditions argument in equilibrium allow for either {statevars + v.X} or {statevars + v.Composition} to be specified? What about {statevars + v.W} (i.e. convert to v.Composition and then to v.X internally)
@bocklund bocklund requested a review from richardotis Dec 11, 2019
@igorjrd

This comment has been minimized.

Copy link
Contributor

igorjrd commented Dec 17, 2019

I think this improvements would be valuable additions to pycalphad. In metallurgy, for example, its common to specify chemical composition in mass fraction instead of mole fraction. Personally, I needed to convert from mass fraction to mole fraction all the times that I used pycalphad to calculate phase diagrams. Also, mix() method is useful to estimate equilibrium or properties in dillution zones and related cases.

At my conversion from mass to mole fraction scripts, I was using a periodic table json database, that I found in pymatgen repository. Maybe this database, or a simplified version of this, can be used to build the masses dictionary. I can contribute with this, if you want some help.

@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Dec 18, 2019

@igorjrd thanks for your feedback, it's really helpful to have input from users!

Regarding the pymatgen database, my main concern if we reuse that file is that the pure element masses that you find in the SGTE unary database are all rounded to 5 significant figures and in some places they even disagree.

For example, Fe in the SGTE unary database is 55.847, but in pymatgen it's 55.845 (which is the accepted value AFAIK). The difference is relatively minor, but by our observation, pycalphad typically agrees with the commercial softwares out to 5 or 6 digits, so this kind of difference can impact whether you see the same answer if you input the same database and conditions into two different software packages.

We have made similar choices in the past to go with a less widely accepted physical constant in choosing the ideal gas constant to be 8.3145 (in agreement with Thermo-Calc), even though the actual value is known with more precision.

@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Dec 23, 2019

This is an excellent design. I only have a few items.

  1. Compositions of species make sense in the context of single phases. Ideally we would support it, but since we don't currently support conditions on the composition of a particular phase, I don't think it's urgent or blocking for this feature.
  2. Dependence on a Database object should be avoided on Composition object creation, if possible. Can conversion to/from mole fractions be made lazy somehow, and can we capture the database context later?
  3. Yes and yes. Probably unpack_conditions should take care of getting arguments into a canonical form.
@bocklund bocklund force-pushed the composition-variables branch from 56d2821 to 654ed4f Dec 23, 2019
@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Dec 30, 2019

@richardotis

  1. The v.Composition objects have optional mass dictionaries now and the mass_fractions/mole_fractions raise errors and prompt the user to set the masses now (either by a dictionary or by hand).

  2. In terms of having an unpack_conditions-like functionality, one issue you could run into is that v.Composition objects as implemented are not compatible with the arange tuples (v.X('A'): (0, 1, 0.01)) or lists/arrays of compositions as we usually unpack. v.Composition objects are strictly point compositions, so mixing them with the ways composition ranges can be specified with unpack_conditions might be odd.

    As the v.Composition objects are written (and since equilibrium doesn't support disabling broadcasting), v.Composition objects would really only be useful when specifying fixed composition, but varying potentials (primarily P, T) equilibrium calculations. If we allow for compositions to not be broadcast (or at least not broadcast for compositions), then v.Compositions become really useful to construct isopleths in compositions.

    My outstanding questions are:

    • What (if any) relationship should v.Composition objects have to ranges of compositions (unpack_conditions style)?
    • Based on the answer to the previous question, should implicit handling and conversion of all v.X/v.W/v.Composition conditions in equilibrium to v.Composition objects be delayed until that feature is more useful (e.g. conditions broadcasting (#172) and/or isopleth plotting (#5))?
@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Dec 30, 2019

  1. ✔️
  2. ✔️
  3. I understand. The main value I see for using a v.Composition object right now would be (a) conversion between mass and mole fraction and (b) isopleth calculation, as you mention. While it's not a blocker for the initial version of the feature, before modifying any of the internals of equilibrium to make implicit Composition conversion work for point compositions, I would like to see the following API made functional (which, I agree, probably involves addressing #172 ):
[...]
comps = ss316.components + ti64.components
conds = {v.P: 101325, v.T: (300, 2000, 20), v.N: 1, **ss316.interpolate_path(ti64, 100)}
eq = equilibrium(dbf, comps, phases, conds)
@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Dec 30, 2019

Also, unrelated to the previous comment and also not a blocker, but I can think of situations where I'd also like to be able to interpolate logarithmic compositions (e.g., interstitial content). Not sure on the right API.

@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Feb 4, 2020

Based on some offline discussion, it was decided to table v.Composition objects since the API didn't seem polished enough.

Instead, this PR will still adopt the v.MassFraction = v.W and v.MoleFraction = v.X as outlined above and two functions in pycalphad.variables exist to convert between them.

Since these conversion functions only work for point compositions, it seems easy enough to convert between point compositions using the v.get_mass_fractions and v.get_mole_fractions functions for use in equilibrium and other code. As it stands, I don't think this has much utility for being worked into equilibrium directly, since it would create so many special cases that would have to be handled, for example if a user tries

  1. broadcasting any v.W (unpack_conditions style) (since the conversions are not broadcastable)
  2. to use v.W conditions mixed with v.X or v.MU conditions
  3. to use v.W('A') and v.X('A') or v.MU('A')

Since complicated composition paths are not very easy to do right now, I think it makes sense to provide the basic functionality for users here and then expand where necessary

@bocklund bocklund marked this pull request as ready for review Feb 4, 2020
Copy link
Collaborator

richardotis left a comment

I really like this API. Great work.

@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Feb 4, 2020

Not sure if this would be considered in scope, but it might be useful for a user who's trying to do calculations with species to be able to convert mole fractions of species to mole fractions of pure elements:

import numpy as np
def get_pure_element_amounts(fractions_dict, dependent_species):
    dependent_species = v.Species(dependent_species)

    indep_pure_element_amounts = {}
    for mole_fraction, amount in fractions_dict.items():
        for pure_element, pure_element_natoms in mole_fraction.species.constituents.items():
            new_amount = (amount*pure_element_natoms/mole_fraction.species.number_of_atoms)
            existing_amount = indep_pure_element_amounts.get(pure_element, 0.0)
            indep_pure_element_amounts[pure_element] = existing_amount + new_amount

    dep_pure_element_frac = 1 - sum(indep_pure_element_amounts.values())
    dep_pure_element_amounts = {}
    for pure_element, pure_element_natoms in dependent_species.constituents.items():
        new_amount = (dep_pure_element_frac*pure_element_natoms/dependent_species.number_of_atoms)
        dep_pure_element_amounts[pure_element] = new_amount
    assert np.isclose(sum(dep_pure_element_amounts.values()), dep_pure_element_frac)

    total_pure_element_amounts = {}
    for pure_element in (indep_pure_element_amounts.keys() | dep_pure_element_amounts.keys()):
        total_pure_element_amounts[pure_element] = indep_pure_element_amounts.get(pure_element, 0.0) + dep_pure_element_amounts.get(pure_element, 0.0)
    assert np.isclose(sum(total_pure_element_amounts.values()), 1.0)
    return total_pure_element_amounts
@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Feb 4, 2020

Can you explain that use case some?

@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Feb 4, 2020

The current PR will let you work in species, but equilibrium only lets you define pure element mole fractions as conditions. The above code will convert v.MoleFractions that are any combination of pure elements and species into the pure element amounts

@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Feb 4, 2020

I don't think in the models that use species that it is generally true that X(A0.3B0.7) = 0.5 == X(A) = 0.5*0.3; X(B) = 0.5*0.7

@richardotis

This comment has been minimized.

Copy link
Collaborator

richardotis commented Feb 4, 2020

Check my math, but when I tried it in Thermo-Calc just now, these two sets of conditions gave different equilibria:

X(AL2O3)=0.1, X(NI)=0.7, N=1, P=1E5, T=2000
X(AL)=4E-2, X(NI)=0.7, T=2000, P=1E5, N=1
@bocklund

This comment has been minimized.

Copy link
Collaborator Author

bocklund commented Feb 4, 2020

X(AL2O3)=0.1, X(NI)=0.7, N=1, P=1E5, T=2000

I can verify the TC stuff as well. I don't understand how this cannot mean there are 0.1*2/5 moles of Al in the system, if there is 0.1 moles of AL2O3 in the system and the dependent component is oxygen

@bocklund bocklund merged commit 2929ba8 into develop Feb 6, 2020
4 of 5 checks passed
4 of 5 checks passed
coverage/coveralls Coverage decreased (-0.08%) to 86.005%
Details
continuous-integration/appveyor/branch AppVeyor build succeeded
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
continuous-integration/travis-ci/push The Travis CI build passed
Details
@bocklund bocklund added this to the 0.8.2 milestone Feb 21, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

3 participants
You can’t perform that action at this time.