Giter VIP home page Giter VIP logo

pycalphad's Introduction

pycalphad, a library for the CALculation of PHAse Diagrams

Join the chat at https://gitter.im/pycalphad/pycalphad

Test Coverage

Build Status

Development Status

Latest version

Supported Python versions

License

Note: Unsolicited pull requests are _happily accepted!

pycalphad is a free and open-source Python library for designing thermodynamic models, calculating phase diagrams and investigating phase equilibria within the CALPHAD method. It provides routines for reading Thermo-Calc TDB files and for solving the multi-component, multi-phase Gibbs energy minimization problem.

The purpose of this project is to provide any interested people the ability to tinker with and improve the nuts and bolts of CALPHAD modeling without having to be a computer scientist or expert programmer.

For assistance in setting up your Python environment and/or collaboration opportunities, please contact the author by e-mail or using the issue tracker on GitHub.

pycalphad is licensed under the MIT License. See LICENSE.txt for details.

Required Dependencies:

  • Python 3.7+
  • matplotlib, numpy, scipy, symengine, xarray, pyparsing, tinydb

Installation

See Installation Instructions.

Examples

Jupyter notebooks with examples are available on NBViewer and pycalphad.org.

Documentation

See the documentation on pycalphad.org.

Getting Help

Questions about installing and using pycalphad can be addressed in the pycalphad Google Group. Technical issues and bugs should be reported on on GitHub. A public chat channel is available on Gitter.

Citing

If you use pycalphad in your research, please consider citing the following work:

Otis, R. & Liu, Z.-K., (2017). pycalphad: CALPHAD-based Computational Thermodynamics in Python. Journal of Open Research Software. 5(1), p.1. DOI: http://doi.org/10.5334/jors.140

Acknowledgements

Development has been made possible in part through NASA Space Technology Research Fellowship (NSTRF) grant NNX14AL43H, and is supervised by Prof. Zi-Kui Liu in the Department of Materials Science and Engineering at the Pennsylvania State University. We would also like to acknowledge technical assistance on array computations from Denis Lisov.

pycalphad's People

Contributors

blazejgrabowski avatar bocklund avatar dependabot[bot] avatar djlicia avatar dschwen avatar gabrielbusta avatar gitter-badger avatar igorjrd avatar jan-janssen avatar jorgepazsoldanpalma avatar jwsiegel2510 avatar kadkison avatar mfrichtl avatar mxf469 avatar nandiniraja348 avatar nury12n avatar richardotis avatar tkphd avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

pycalphad's Issues

Muggianu correction is incorrectly applied in higher-order systems

It looks like the Muggianu ternary excess correction routine is called inside of a loop, so it is recursively substituting on itself. This problem only shows up in systems with more than 3 components which have ternary interaction parameters.

The fix is to move the correction to be applied to a local variable.

problem with calculating equilibirum

Hi Otis,

Recently, I tried to do some calculation which relies on the computation of equilibrium for the Al-Ni system with FCC_A1 and gamma prime phase. I tried to do this with the pycalphad, but it turned out that there are some difference among different implementations of the computation of equilibrium.

By setting temperature as 873.15K, pressure as 101325, and the molar fraction of Al as 0.2, I got results from Thermo-Calc, opencalphad and pycalphad.

Result from opencalphad(the result from Thermo-Calc is the same as opencalphad.)


Some data for components ...................................:
Component name    Moles      Mole-fr  Chem.pot/RT  Activities  Ref.state
AL                2.0000E-01  0.20000 -2.4583E+01  2.1080E-11  SER (default)   
NI                8.0000E-01  0.80000 -5.2649E+00  5.1701E-03  SER (default)   

Some data for phases .......................................:
Name                Status Moles      Volume    Form.Units  At/FU dGm/RT  Comp:
FCC_A1.................. E  2.626E-01  0.00E+00  2.63E-01    1.00  0.00E+00  X:
 NI     8.94297E-01  AL     1.05703E-01
 Constitution: Sublattice  1 with     2 constituents and     1.000000 sites
 NI     8.94297E-01  AL     1.05703E-01
               Sublattice  2 with     1 constituents and     1.000000 sites
 VA     1.00000E+00

GAMMA_PRIME............. E  7.374E-01  0.00E+00  7.37E-01    1.00  0.00E+00  X:
 NI     7.66422E-01  AL     2.33578E-01
 Constitution: Sublattice  1 with     2 constituents and     0.750000 sites
 NI     9.99841E-01  AL     1.58647E-04
               Sublattice  2 with     2 constituents and     0.250000 sites
 AL     9.33838E-01  NI     6.61624E-02

Result from pycalphad


X = [ 0.2272507,  0.7727493], [ 0.1060816,  0.8939184]
Y = [  3.64443392e-04,   9.99635557e-01,   9.07909485e-01, 9.20905148e-02],
          [ 1.06081603e-01,   8.93918397e-01,   1.00000000e+00, nan]

At present, I would like to use opencalphad because it produces the same result as the Thermo-Calc, though the result from pycalphad is closed to the opencalphad. I wonder whether this is problem for pycalphad. I report this and hope this might be helpful.

By the way, the APIs designed in pycalphad is very convenient. And here is how I called pycalphad to obtain the result above


db = Database('alnicrre.tdb')
my_phases = ['FCC_A1', 'GAMMA_PRIME']
eq = equilibrium(db, ['AL', 'NI', 'VA'], my_phases, {v.X('AL'): 0.2, v.T: 873.15, v.P: 101325})
print eq.get("X")
print eq.get("Y")

Best wishes,

Jing

Boolean indexing errors while computing equilibrium

The symptom is VisibleDeprecationWarnings about boolean indexing. The cause is the simplex filter in lower_convex_hull doesn't properly handle the case of when two or more different trial simplices for the same condition satisfy the non-negativity criteria, i.e., literally an edge/corner case. I've only seen it when you try to compute right on top of stoichiometry for a B2 or L12 phase. In principle, even if you get a result the answer is not guaranteed to be correct because the arrays will have been misaligned. I think this will happen any time we try to compute equilibrium right on top of a sublattice end-member.

dbf_dupin = Database('Al-Ni/Al-Ni-Dupin-2001.tdb')
ni3al_dupin = equilibrium(dbf_dupin, ['NI', 'AL', 'VA'], 'BCC_B2', {v.T: (300, 1600, 100), v.P: 101325, v.X('AL'): .5})

You can get around it by changing v.X('AL') to .5001 or similar.

tdb: reading/writing diffusion data

I received a request recently to add support for reading diffusion parameters from TDB files. This looks like it would be straightforward and would only require a slight modification to the parser for the PARAMETER keyword to allow the &component syntax for diffusing species. Beyond that, the Database object should be flexible enough to allow users to write a Model subclass to construct diffusion coefficients using a custom query passed to the redlich_kister_sum function.

Sorting elements in Redlich-Kister interaction parameters

It seems that most (all?) CALPHAD software out there is alphabetically sorting the components listed in binary (and higher-order?) interaction parameters. What this means is that, e.g., G(LIQUID,C,W;1) == G(LIQUID,W,C;1) even though they theoretically differ by a sign when the order parameter is odd (1 in this case). I haven't been able to find this documented anywhere but it has an impact on pycalphad's ability to reproduce results from the literature, e.g. the W-C diagram from Gustafson 1986, and I'm sure there are many others.

Installation and examples

Dear Richard,

The pycalphad is really great project, so I wanted to test it out. I do have problems with the pycalphad examples. I am running Mac Pro with Python 2.7.8 (default, Nov 12 2014, 14:36:21)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.54)] on darwin, with standard scipy stack.
When I run BinTest for example I get:

ImportError                               Traceback (most recent call last)
<ipython-input-9-c9b5ef89e019> in <module>()
      1 get_ipython().magic(u'matplotlib inline')
----> 2 from pycalphad import Database, binplot
      3 
      4 db = Database('alzn_mey.tdb')
      5 my_phases = ['FCC_A1', 'HCP_A3', 'LIQUID']

ImportError: cannot import name binplot

It seems that binplot is not found. I did installation with: python setup.py install from current github repository.

Thank you very much for your help, and keep up the great work!
Best regards,
โˆƒt

Order-disorder transition lines

Hi Richard, it would be good to implement a order-disorder transition plotting tool that would look at the differences in the site fractions. For example in a 2 sublattice BCC_B2, model, it would calculate when x(comp1) in site 1 is different from x(comp 1) in site 2. This has to be done manually in Thermo-Calc right now and sometimes it does not work correctly without trying many compositions first.

Chemical potential, driving force and changing model-parameters

Hello Richard,
I would like to use pycalphad for a simulation.
Each step requires:

  1. calculation of chemical potential and driving forces for phases.
  2. changing parameters for the Gibbs energy description of some phases.
  3. making equilibrium calculation

Now, for (1), is there a simple means of calculating chemical potential of a component in phase ? Can you give an example?

For (2) I suspect that the answer could come from the fitting module. If I include some adjustable parameters in the description of a phase, changing these parameters may be a way to modify its Gibbs energy. Can this be done in the current version of pycalphad?

Thanks,
Eli Brosh

Mixing properties are incorrect for phases with non-chemical energy contributions

The reason is, e.g., GM_MIX, is only subtracting out the contribution from self.models['ref'], which only includes the chemical energy contribution to each end-member. For phases with magnetic, Einstein, or some custom end-member contributions, the result will not be correct.

Standard properties like GM, CPM, etc., appear to be fine.

Order-disorder model fails to build when VA is not an active component

This is the unit test failure:

check_energy(Model(DBF, ['AL'], 'B2'), {v.T: 1400, v.SiteFraction('B2', 0, 'AL'): 1,
v.SiteFraction('B2', 1, 'AL'): 1}, -6.57639e4, mode='sympy')
======================================================================
ERROR: Pure component end-members in sympy mode.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/rotis/.virtualenvs/calphad/lib/python3.3/site-packages/nose/case.py", line 198, in runTest
    self.test(*self.arg)
  File "/home/rotis/git/pycalphad/pycalphad/tests/test_energy.py", line 298, in test_pure_sympy
    check_energy(Model(DBF, ['AL'], 'B2'), \
  File "/home/rotis/git/pycalphad/pycalphad/model.py", line 44, in __init__
    phase.upper(), dbe.symbols, dbe.search)
  File "/home/rotis/git/pycalphad/pycalphad/model.py", line 147, in build_phase
    symbols, param_search)
  File "/home/rotis/git/pycalphad/pycalphad/model.py", line 586, in atomic_ordering_energy
    molefraction_dict[sitefrac] = species_dict[sitefrac.species]
KeyError: 'VA'

This doesn't happen in the other models because they all check their contributions against Model.components to ensure that only active components are included.

Plotting arbitrary isopleths

Users would have a plot() style interface for specifying arbitrary sections through composition space. This is not too conceptually difficult using the energy_surf routine and scipy's ConvexHull routine (though we may run into performance problems); the main trick is efficiently calculating intersections of hyperplanes with the convex hull. Routines exist for this but need to be investigated further.

Apparent error in calculation of driving force

Hi Richard,
I installed pycalphad 0.2.1
pycalphad.version
'0.2.1'

Now I tried to calculate the driving force for precipitation of the Al13Fe4 phase from liquid Al-Fe at x(Fe)=0.1

Here is the script:

from pycalphad import Database, binplot, equilibrium
import pycalphad.variables as v
from numpy import *
from pylab import *

db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')

R=8.31451
T=arange(1000,2001,20)
xAl=.9
xFe=1-xAl
data = equilibrium(db, ['AL', 'FE'], 'LIQUID', {v.X('AL'): xAl, v.T: T, v.P: 1e5},verbose=False)

muFe=data['MU'].sel( component='FE').data[0].flatten()
muAl=data['MU'].sel( component='AL').data[0].flatten()

data = equilibrium(db, ['AL', 'FE','VA'], 'AL13FE4', {v.X('AL'): 13/17.0, v.T: T, v.P: 1e5},verbose=False)
G=data['GM'].data[0].flatten()

dGM=-(G-xAl_muAl-xFe_muFe)/R/T

figure(1)
plot(T,dGM,'g')
dgm

show(block=False)

The calculation gives positive driving force below 1800 K.
It should be positive below ~1250 K as can be seen on the Al-Fe phase diagram.
A corresponding result of a thermo-calc calculation is also attached (with the correct result).
dgm_thermocalc

It seems that the calculation of the Gibbs energy of Al13Fe4 is erroneous in pycalphad

For reference, Gm(Al13Fe4) mu(Al) and mu(Fe) at T=2000 x(Al)=0.9 x(Fe)=.1 are given

pycalphad:
mu(Al)=-130684.44533461
mu(Fe)=-202444.57089581
Gm(Al13Fe4)=-134716.70446316

Thermo-calc:
mu(Al)=-130684.45
mu(Fe)=-202444.6
Gm(Al13Fe4)=-136221.09

Eli

equilibrium: Detecting single-phase regions

Hi Richard,
I get a strange result from a seemingly simple calculation:

from pycalphad import Database, equilibrium
import pycalphad.variables as v

db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')

data = equilibrium(db, ['AL', 'FE','VA'], ['LIQUID','AL13FE4'], {v.X('FE'): 1.0e-1, v.T: 1500, v.P: 1e5},verbose=False)

print('phases')
print(data['Phase'].values.flatten())
print(' ')
print('phase fractions')
print(data['NP'].values.flatten())
print(' ')
print('composition of the two phases')
print(data['X'].sel(component='FE').values.flatten())

The output is:
phases
['LIQUID' 'LIQUID']

phase fractions
[ 0.58441188 0.41558812]

composition of the two phases
[ 0.09999928 0.10000101]

Why two liquid phases instead of one?
There should not be a miscibility gap there.

Best regards,
Eli

binplot sometimes hangs on Windows

I'm pretty sure this is caused by the underlying qhull library (provided by scipy) not handling large arrays of points very well. I'm not sure why Linux/OSX seems to behave better. It's expected that once the routine is updated to use the new equilibrium engine (the same one used by equilibrium), this bug will go away, since qhull is not used there.

Proposed syntax for equilibrium calculation

These are some sketches of ideas for a unified user interface for all the equilibrium computation in pycalphad. The goal is to implement this interface for the next release (0.2).

import numpy as np
import pycalphad.variables as v
from pycalphad import Database, equilibrium

Point calculation

dbf = Database('database.tdb')
phases = ['PHASE1', 'PHASE2', 'PHASE3']
equilibrium(dbf, ['A', 'B', 'VA'], phases, {v.X('A'): 0.55, v.N: 1, v.T: 1400},
    pdens=1000)

Note: Need to add routine for unpacking conditions, including

  • Convert strings to pycalphad.variables objects
  • Convert conditions to starting values for all phases

Setting phase status before a calculation

phases = {'PHASE1': 'suspended', 'PHASE2': 'dormant', 'PHASE3': 'entered'}

Suspended: Not considered in the calculation at all (default if not included in a list)

Dormant: Included for the purposes of driving force computation, but not allowed in any equilibria.
The way to implement this is to perform a separate lower convex hull computation for this phase along the calculation path.

Entered: Fully considered in the calculation (default if included in a list)

Step calculation

equilibrium(dbf, ['A', 'B', 'VA'], phases, {v.X('A'): 0.55, v.N: 1, v.T: (300,)})
equilibrium(dbf, ['A', 'B', 'VA'], phases, {v.X('A'): 0.55, v.N: 1, v.T: (300, 1400)})
equilibrium(dbf, ['A', 'B', 'VA'], phases, {v.X('A'): 0.55, v.N: 1, v.T: (300, 1400, 5)})
equilibrium(dbf, ['A', 'B', 'VA'], phases,
                 {v.X('A'): 0.55, v.N: 1, v.T: np.linspace(300, 1400, num=100)})

Rules for keys of conditions dicts:

  • If it's numeric, treat as a point value
  • If it's a tuple with one element, treat as a point value
  • If it's a tuple with two elements, treat as lower/upper limits and guess a step size
  • If it's a tuple with three elements, treat as lower/upper/step
  • If it's a list, ndarray or other non-tuple ordered iterable, use those values directly

Map calculation

equilibrium(dbf, ['A', 'B', 'VA'], phases,
                 {v.X('A'): (0, 1, 0.01), v.N: 1, v.T: (300, 1400, 5)})

equilibrium(dbf, ['A', 'B', 'C', 'VA'], phases, {v.X('A'): (0, 1, 0.01),
                                                 v.X('B'): (0, 1, 0.01),
                                                 v.N: 1, v.T: (300, 1400, 5)})

equilibrium(dbf, ['A', 'B', 'C', 'VA'], phases, {v.X('A'): (0, 1, 0.01),
                                                 v.X('C'): np.logspace(-12, -1, num=20),
                                                 v.N: 1, v.T: (300, 1400, 5)})

Custom path calculation

The key idea is to accept a list of conditions dicts.

helical_path = [{v.X('A'): np.sin(t),
         v.X('B'): np.cos(t),
         v.N: 1, v.T: (300, 1400)} for t in np.linspace(0, 2*np.pi)]
equilibrium(dbf, ['A', 'B', 'C', 'VA'], phases, helical_path)

Return values

Return an ndarray of EquilibriumResult with the same shape as broadcasting all conditions together.
If a list of conditions dicts is passed, return a list of EquilibriumResults.

If a particular set of conditions is infeasible, determined either by the composition conditions (e.g., fractions sum to more than one) or through the computation (e.g., activity condition cannot be met or system is under/overconstrained),
return an EquilibriumResult with the specified conditions, but with Nones for values.

Note: Given how many of these objects might be made, should we implement EquilibriumResult with __slots__ to try to reduce the memory overhead?

Remark: Would this be a good candidate for the xray project? It needs to be easy to retrieve results for specific conditions, or sets of conditions. Could we get away with just using pandas and storing the results in a DataFrame? We probably need real metrics for this.

How do we postprocess the EquilibriumResults into specific plots (property diagrams, phase diagrams)? One thought is we should leave this to plot-specific routines that will drive the equilibrium computation and handle the details.

Return values (alternative)

What if we got rid of EquilibriumResult for non-point calculations, and returned DataFrames with T, P, energy, entropy, composition, chemical potentials and heat capacities? What about adding an out kwarg to specify which output values we were interested in, so an isothermal phase diagram calculation, where we need energies of each phase and chemical potentials to draw the diagram, might involve:

equilibrium(dbf, ['A', 'B', 'C', 'VA'], phases,
            {v.X('A'): (0, 1), v.X('B'): (0, 1), v.N: 1, v.T: 1200},
            out=[v.GM, v.MU])

This still isn't a fully formed idea, but it seems promising and could allow us to eliminate the tricky function energy_surf from the user-facing API.

Left to a future release

  • Activity conditions
  • Chemical potential conditions
  • Phase amount conditions

Phase diagram plotting generates spurious tie lines

This bug occurs when the phase diagram plotter has trouble distinguishing a true tie line resulting from a miscibility gap from two points in a single-phase region. The cause of this bug is two things: (1) Not enough points sampled in that region and/or (2) The minimum tie line length is set too small.

Fixing this bug requires rewriting the point sampling routine to guarantee a minimum point density to the plotter everywhere in the diagram. The plotter routines will then know that "tie lines" shorter than the minimum nearest-neighbor distance are not real.

Example:
image

Platform-dependent binplot crash bug

I'm pretty sure this is the same bug Jeff reported by e-mail a few weeks ago, which means the solution of filtering out nan's will probably work as a stopgap before a total rewrite.

ReadTheDocs is failing to build

https://readthedocs.org/projects/pycalphad/builds/3465916/

Installed /home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/lib/python2.7/site-packages/tinydb-2.4-py2.7.egg
Searching for autograd
Reading https://pypi.python.org/simple/autograd/
Best match: autograd 1.1.2
Downloading https://pypi.python.org/packages/source/a/autograd/autograd-1.1.2.tar.gz#md5=db1cc306f37bb39f30d36c72a39ce1e9
Processing autograd-1.1.2.tar.gz
Writing /tmp/easy_install-Wtdg_w/autograd-1.1.2/setup.cfg
Running autograd-1.1.2/setup.py -q bdist_egg --dist-dir /tmp/easy_install-Wtdg_w/autograd-1.1.2/egg-dist-tmp-Sj6EJJ
Traceback (most recent call last):
  File "setup.py", line 47, in <module>
    'Programming Language :: Python :: 3.4'
  File "/usr/lib/python2.7/distutils/core.py", line 151, in setup
    dist.run_commands()
  File "/usr/lib/python2.7/distutils/dist.py", line 953, in run_commands
    self.run_command(cmd)
  File "/usr/lib/python2.7/distutils/dist.py", line 972, in run_command
    cmd_obj.run()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/install.py", line 67, in run
    self.do_egg_install()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/install.py", line 117, in do_egg_install
    cmd.run()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 380, in run
    self.easy_install(spec, not self.no_deps)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 610, in easy_install
    return self.install_item(None, spec, tmpdir, deps, True)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 661, in install_item
    self.process_distribution(spec, dist, deps)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 709, in process_distribution
    [requirement], self.local_index, self.easy_install
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 836, in resolve
    dist = best[req.key] = env.best_match(req, ws, installer)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 1081, in best_match
    return self.obtain(req, installer)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 1093, in obtain
    return installer(requirement)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 629, in easy_install
    return self.install_item(spec, dist.location, tmpdir, deps)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 659, in install_item
    dists = self.install_eggs(spec, download, tmpdir)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 842, in install_eggs
    return self.build_and_install(setup_script, setup_base)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 1070, in build_and_install
    self.run_setup(setup_script, setup_base, args)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/easy_install.py", line 1056, in run_setup
    run_setup(setup_script, args)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 240, in run_setup
    raise
  File "/usr/lib/python2.7/contextlib.py", line 35, in __exit__
    self.gen.throw(type, value, traceback)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 193, in setup_context
    yield
  File "/usr/lib/python2.7/contextlib.py", line 35, in __exit__
    self.gen.throw(type, value, traceback)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 164, in save_modules
    saved_exc.resume()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 139, in resume
    compat.reraise(type, exc, self._tb)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 152, in save_modules
    yield saved
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 193, in setup_context
    yield
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 237, in run_setup
    DirectorySandbox(setup_dir).run(runner)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 267, in run
    return func()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 236, in runner
    _execfile(setup_script, ns)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/sandbox.py", line 46, in _execfile
    exec(code, globals, locals)
  File "/tmp/easy_install-Wtdg_w/autograd-1.1.2/setup.py", line 79, in <module>

  File "/usr/lib/python2.7/distutils/core.py", line 151, in setup
    dist.run_commands()
  File "/usr/lib/python2.7/distutils/dist.py", line 953, in run_commands
    self.run_command(cmd)
  File "/usr/lib/python2.7/distutils/dist.py", line 972, in run_command
    cmd_obj.run()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/bdist_egg.py", line 151, in run
    self.run_command("egg_info")
  File "/usr/lib/python2.7/distutils/cmd.py", line 326, in run_command
    self.distribution.run_command(command)
  File "/usr/lib/python2.7/distutils/dist.py", line 972, in run_command
    cmd_obj.run()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/egg_info.py", line 180, in run
    self.find_sources()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/egg_info.py", line 207, in find_sources
    mm.run()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/egg_info.py", line 291, in run
    self.add_defaults()
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/egg_info.py", line 320, in add_defaults
    sdist.add_defaults(self)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pycalphad/envs/develop/local/lib/python2.7/site-packages/setuptools/command/sdist.py", line 130, in add_defaults
    build_ext = self.get_finalized_command('build_ext')
  File "/usr/lib/python2.7/distutils/cmd.py", line 312, in get_finalized_command
    cmd_obj.ensure_finalized()
  File "/usr/lib/python2.7/distutils/cmd.py", line 109, in ensure_finalized
    self.finalize_options()
  File "/tmp/easy_install-Wtdg_w/autograd-1.1.2/setup.py", line 22, in finalize_options
    url='https://github.com/richardotis/pycalphad',
AttributeError: 'dict' object has no attribute '__NUMPY_SETUP__'

limiting the composition range of a phase

Hello Richard,
I am looking for a way to make a calculation with a constraint that a certain phase can appear only at a defined composition range.
Suppose we have two phases with Gibs energy -composition curves as depicted in the attached figure.
energy curves

Suppose the overall composition is x0=0.65
An equilibrium calculation for this system would give a two-phase coexistence.
I would like to be able to constrain the calculation so that phase two cannot appear at certain compositions. For example, If I would be able to say that phase 2 cannot appear with composition x>0.6, the constrained equilibrium calculation would give a single phase equilibrium (only phase 1).
This is equivalent to adding a very large positive value to the Gibbs energy of phase 2 for x>0.6:
G2=G2+(x>0.6)*1e10.

Is there a way to do this in pycalphad?

Regards,
Eli

xray dependency renamed to xarray

See pydata/xarray#704.

According to their plan the xray and xarray namespaces will coexist for a little while, but the former will be deprecated as of v0.7. Once that release is out, we will push the rename.to users on our end.

  • xarray 0.7 released
  • Rename xray to xarray in all imports
  • Update setup.py
  • Update conda recipe
  • Update unit tests
  • Update CI configuration files

Numerical stability issue with equilibrium calculation

Hello Richard,
Shana Tova.
I tried to calculate entropy of formation of a compound as function of composition.
The calculation of entropy is done by numerical differentiation of Gibbs energy of formation.
See code below.
The result, as can be seen in the attached figure is noisy.
Perhaps this can be corrected by calculating the entropy not by numerical differentiation.
Is there a way to use symbolic algebra for this?

Thanks,
Eli

from pycalphad import Database, equilibrium
import pycalphad.variables as v
from numpy import *
from pylab import *




db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')



def sb(T,phase_name,x_liquid,x_phase,db):
    dT=1.0

    f1m=dG(T-1*dT,phase_name,x_liquid,x_phase,db)
    f0=dG(T+0*dT,phase_name,x_liquid,x_phase,db)

    DG=f0-f1m
    s=DG/(1.0*dT)
    print(T,x_liquid,s)
    return s



def dG(T,phase_name,x_liquid,x_phase,db):
        data = equilibrium(db, ['AL', 'FE'], 'LIQUID', {v.X('FE'): x_liquid, v.T: T, v.P: 1e5},verbose=False)
        muFe=data['MU'].sel( component='FE').data[0].flatten()
        muAl=data['MU'].sel( component='AL').data[0].flatten()
        data = equilibrium(db, ['AL', 'FE','VA'], phase_name, {v.X('FE'): x_phase, v.T: T, v.P: 1e5},verbose=False)
        G=data['GM'].data[0].flatten()
        dG=(G-(1-x_phase)*muAl-(x_phase)*muFe)

        return dG


n=100.0
irange=arange(0,n+1,1)    
xrange=irange/n
xrange[0]=1.0e-4
xrange[-1]=.9999
T=1000*ones(size(irange))
S=zeros(size(irange))

for i in irange: 
    S[i]=sb(T[i],'AL13FE4',xrange[i],4/17.0,db)


plot(xrange,S)
xlabel('x(liquid,Fe)')
ylabel('Delta(S) formation of AL13FE4' )


show(block=False)

noise in entropy of formation

Deprecating Python 3.3 support

Upstream is no longer building conda packages for Python 3.3, and that is making it difficult to retain and test support for that Python version. The xarray package rename is really what forces the issue, since changing the dependency resolution has exposed how far behind py33 really is.

We probably need to discontinue support soon.

Drop Python 3.3 support

It looks like building "clean" conda packages supporting Python 3.3 isn't possible anymore without leaning on pip, which is what we do on Travis CI. We can still test against it, but it's probably best to say 3.4 is the minimum supported version now.

Ternary isothermal phase diagram plotting

Ternary isothermal diagrams have a lot of performance and accuracy problems, especially when dealing with phases with internal degrees of freedom but small solubility. See the below example for Cr-Ti-V.

image

Correct diagram (from NIMS phase diagram database):
image

The complete solution isn't clear, but it probably starts with a binplot-style rewrite and heavy CPU/memory profiling. Verification of the test databases using commercial software should also be done prior to any further work.

Release on conda-forge

It looks pretty easy to add recipes to conda-forge, so this should probably be the preferred way to release pycalphad using conda in the future. For the 0.4 release we will simultaneously release on PyPI and the richardotis conda channel as per usual, but we will also submit a new recipe for pycalphad 0.4 to conda-forge. If the process goes well, we will deprecate the richardotis channel and do all future conda releases exclusively on conda-forge.

Getting output of Gibbs energy, chemical potential, atomic mobility, etc.

From @lightinight

Hello Otis,

The pycalphad is a very interesting project. As a matter of fact, I'm writing a similiar parser to read a tdb file and output the expressions for molar Gibbs energy, chemical potential, thermodynamic factor, atomic mobility and so on. They are very essential to phase-field modeling and simulation. I wonder whether or not you are interested in such work. And I would be very appreciated if I could some contribution to such work.
At present, I have implemented a very basic version to provide those thermodynamic and kinetic infomation for phases with a single sublattice. I have no idea how the Gibbs energy of a multi-sublattice phase should be constructed.

Best wishes!

Zhong JIng

TDB parser doesn't understand % symbol in component names

The parser for TDB files does not properly handle the "%" symbol in component names. This causes it to silently fail to parse CONSTITUENT commands containing these symbols, causing opaque error messages for the end user.

Workaround

Delete % symbols from component names.

Solution

The parser should drop these symbols silently since pycalphad's equilibrium engine won't need hints on the starting point for energy minimization.

Invalid number of arguments error when energy simplifies to zero

It's possible to get a difficult-to-debug error concerning the number of arguments passed to the energy function. This happens if the energetic contribution simplifies to zero. For example:

PHASE M7C3_D101 %  2 7 3 !
CONSTITUENT M7C3_D101 :MN:C: !
PARAMETER G(M7C3_D101,MN:C;0)  1 VV22+VV23*T**2+VV24*T**3
 +VV25*T**4+VV26*T**5; 6000 N !

where all of the above VV variables are zero. If the user happens to typo the name of the phase in the TDB, they may be scratching their head about this error (perhaps parameters without a matching phase should throw their own error, but that's a separate issue).

Chemical potentials calculated by pycalphad are discontinuous

By means of pycalphad, the curved surface of the chemical potentials calculated are discontinuous. As is shown in the figure1, this is the scatter diagram of the chemical potentials of component Ni of the Al-Ni-Cr ternary system at 1273K. In the figure1, some of the chemical potentials change dramatically in different compositions.
figure1
Selecting a composition corresponding to the dramatically variational chemical potential from figure1, calculate the chemical potentials of different compositions around that composition. As is shown figure2, this is the diagram of the chemical potentials of component Ni of compositions of the area.
figure2
How to solve the problem that the chemical potentials are discontinuous calculated by pycalphad?

Chemical potentials calculated by pycalphad are dis.docx

"zip_longest" vs "izip_longest" issue

I had to modify "itertools.zip_longest" on line 144 of "pycalphad-master\pycalphad\plot\binary.py" to "itertools.izip_longest" to work with python 2.7.8.

From the itertools documentation, it looks like maybe this is a difference between the python 2 and python 3 versions of itertools.

equilibrium: CalledProcessError on Windows

Many users on Windows do not have the MSVC C/C++ compiler installed, causing the compiled backend to break with cryptic error messages. It does not appear to be easy to distribute or otherwise provide users with this compiler. The best option may be to do some auto-detection or error checking for the build process, and then catch CalledProcessError and fall back to the interpreted backend (with a warning). We can provide instructions in the documentation for how to improve performance on Windows by installing MSVC, either the Python 2.7 compiler edition or the MSVC 2015 Community Edition on Python 3.5.

Cleaning up the API for custom Models

To add custom energetic contributions to the default Model class, you currently have to do something like the following:

from pycalphad import Model
from tinydb import where

# Necessary for TDB compatibility
import pycalphad.io.tdb_keywords
pycalphad.io.tdb_keywords.append('CUSTOM')

class CustomModel(Model):
    def custom_energy(phase, param_search):
        param_query = (
            (where('phase_name') == phase.name) & \
            (where('parameter_type') == "CUSTOM") & \
            (where('constituent_array').test(self._array_validity))
        )
        energy_term = self.redlich_kister_sum(phase, param_search, param_query)
        return energy_term / self._site_ratio_normalization(phase)
    def build_phase(self, dbe, phase_name, symbols, param_search):
        """
        Apply phase's model hints to build a master SymPy object.
        """
        phase = dbe.phases[phase_name]
        self.models['ref'] = self.reference_energy(phase, param_search)
        self.models['idmix'] = self.ideal_mixing_energy(phase, param_search)
        self.models['xsmix'] = self.excess_mixing_energy(phase, param_search)
        self.models['mag'] = self.magnetic_energy(phase, param_search)

        # ACTUAL CUSTOM CODE HERE
        self.models['custom'] = self.custom_energy(phase, param_search)
        # END CUSTOM CODE

        # Next, we handle atomic ordering
        # NOTE: We need to add this one last since it uses self.models
        # to figure out the contribution.
        # It will also MODIFY self.models
        ordered_phase_name = None
        disordered_phase_name = None
        try:
            ordered_phase_name = phase.model_hints['ordered_phase']
            disordered_phase_name = phase.model_hints['disordered_phase']
        except KeyError:
            pass
        if ordered_phase_name == phase_name:
            self.models['ord'] = self.atomic_ordering_energy(dbe,
                                                             disordered_phase_name,
                                                             ordered_phase_name)

There is a lot of boilerplate code, and things like writing parameter queries exposes a lot of the internals of the class and requires using Model's internal API. The purpose of this issue is to track brainstorming a cleaner API for

  1. Adding custom energetic contributions while still ensuring correct order so things like order-disorder contributions still work
  2. Adding custom properties, e.g., TC
  3. Writing Database parameter queries (maybe this should be a Database convenience method?)

a switch in chemical potential calculation

I calculated the chemical potential in of Fe in Al-Fe liquid.
Here is the code:
db = Database('alfe_sei.TDB')
data = equilibrium(db, ['FE', 'AL'], 'LIQUID', {v.X('FE'): .2, v.T: Trange, v.P: 1e5},verbose=False)
muFe=data['MU'].sel( component='FE').data[0][0][0]
print('Mu(Fe)=',muFe)

The output is
MuFe=-133484.35563831546

Now, this is not what I got from thermo-calc with the same calculation.
Interestingly, the value of -133484.356 is the chemical potential of Al for the same conditions.
So, perhaps pycalphad somehow switched between the two potentials.

Best regards,
Eli

New project domain: pycalphad.org

https://readthedocs.org/projects/pycalphad/builds/4264743/
https://readthedocs.org/projects/pycalphad/builds/4264744/

Relevant error:

error: python-dateutil 1.5 is installed but python-dateutil<3.0.0,>=2.1 is required by set(['botocore'])

This might be a good excuse to move off of RTD for good, especially since it is still before the first pycalphad paper has been submitted (so we can change the docs URL etc.)

@bocklund Have you made any progress with Travis-CI and gh-pages integration?

the precision of the values of the Gibbs energy

Hi,

Recently, I did some comparison of the values of the Gibbs energy from Tabulate module of the Thermo-Calc, pycalphad and my reader, which shows there are slightly difference between the values.

For the Gibbs energy of the FCC_A1 in the Al-Ni system(with magnetic ordering), the conditions are set as T = 873.15K, p = 101325, x(al) = 0.01 and the values are listed as following

  • TC: -3.84128821E+04 (R = 8.31451)
  • pycalphad: -38412.9016962598 (R = 8.31451) and -38412.9011777456(R = 8.3145)
  • mine: -38412.8820655441 (R = 8.31451)

It seems that the value of the gas constant do make a difference. But such difference might not be as larger as the one between the pycalphad and the TC.

Just report this and hope it might be helpful.

Best wishes!

Release blockers for 0.2

  • Installation instructions for Windows, OSX, and Linux
  • A build for Anaconda so Windows users can conda install pycalphad
  • Simple examples for the documentation, copied from research and example notebooks

Regression on numpy 1.10

It looks like some stricter type casting conversion rules for ufuncs in numpy 1.10 has broken a whole bunch of code which uses np.linalg.solve and np.einsum. Some objects which I thought were of dtype float were apparently actually generic objects getting implicitly converted. It will take some time to fix.

Ship tests/devtools as extras

Per a suggestion from @bocklund it would be good to add a tests extra to make installing the tests easier for benchmarking purposes, so that it's possible to run tests outside of development mode.

One challenge is I don't see an easy way to ship pycalphad.tests only when the extra is specified.

Any solution needs to be added to setup.py as well as our conda recipe (though that process also needs to change now that we're on conda-forge; see #45 )

Regression on pyparsing 2.1.0

Database loading appears to hang on pyparsing 2.1.0.

Reporting user traceback:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pycalphad/io/database.py", line 110, in __new__
    return cls.from_file(fname, fmt=fmt)
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pycalphad/io/database.py", line 183, in from_file
    format_registry[fmt.lower()].read(dbf, fd)
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pycalphad/io/tdb.py", line 747, in read_tdb
    tokens = _tdb_grammar().parseString(command)
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1152, in parseString
    loc, tokens = self._parse( instring, 0 )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2554, in parseImpl
    ret = e._parse( instring, loc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2440, in parseImpl
    loc, exprtokens = e._parse( instring, loc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2695, in parseImpl
    return self.expr._parse( instring, loc, doActions, callPreParse=False )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2423, in parseImpl
    loc, resultlist = self.exprs[0]._parse( instring, loc, doActions, callPreParse=False )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2695, in parseImpl
    return self.expr._parse( instring, loc, doActions, callPreParse=False )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1083, in _parseCache
    value = self._parseNoCache( instring, loc, doActions, callPreParse )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1018, in _parseNoCache
    loc,tokens = self.parseImpl( instring, preloc, doActions )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 2625, in parseImpl
    if e.canParseNext(instring, tmpLoc):
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1066, in canParseNext
    self.tryParse(instring, loc)
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1060, in tryParse
    return self._parse( instring, loc, doActions=False )[0]
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 1080, in _parseCache
    return (value[0],value[1].copy())
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 548, in copy
    ret = ParseResults( self.__toklist )
  File "/Users/firebird/allPythonEnv/pycalphad031/lib/python2.7/site-packages/pyparsing.py", line 285, in __init__
    self.__asList = asList
KeyboardInterrupt

equilibrium: Vectorizing the refinement step

We're leaving probably a 10-100x speedup on the table by not vectorizing the refinement (aka "Solve (3/3)") step in the equilibrium calculation. It's not reasonable to expect users to wait an hour or more to perform mapping calculations or compute a phase diagram. The challenge is that the current serial algorithm is much, much easier to reason about and it will take some careful thinking to get all the indexing operations correct.

Here is a raw dump of notes about how it could work.

First I need to initialize the Hessian matrix and the b vector
Let's try some very conservative upper bounds
max_phases = len([c for c in comps if c != 'VA'])
max_num_vars = max_internal_dof*max_phases+max_phases
max_constraints = num_mass_bals + max_sublattices*max_phases
CONSIDER THAT we need to include an extra 'max_phases' dimension because fancy indexing is done in parallel
A phase participating in the same equilibrium multiple times will index the same condition simultaneously -- this will break
The solution is to use the phase_idx as an extra dimension, then sum over that dimension before copying it over to l_hessian
Hessian size will be (*conds, max_phases, max_num_vars+num_constraints, max_num_vars+num_constraints)
b vector size will be (*conds, max_phases, max_num_vars+num_constraints)
Initialize l_multipliers as zeros of shape (*conds, num_constraints)
Set l_multipliers[*conds, :len(max_phases)] = chemical potentials
Initialize Hessian as stacked identity matrix
Initialize b vector as stacked zero vector
Now we need to fill out the Hessian and b at all conditions, for the current iteration
We know the phases participating at all conditions
We know the best guess site fractions for all participating phases at all conditions
For each phase, make a flattened list of (P, T, *sitefracs) to compute -- filter out converged conditions
   Also track the conditions multi-index for each phase computation
   Also track the var_idx we can write sitefrac-related dof to, and also the phase_idx for phasefrac-related dof
   So it should probably be three flattened lists of equal length
   Pretty easy to compute G, G' and G'' for the flattened list -- where do we copy the results to?
   first_var_idx = index of first site fraction variable for a given phase
   last_var_idx = index of last site fraction variable for a given phase
   Hessian matrix (don't forget to clip P and T derivatives):
   Sitefrac/sitefrac derivatives
     Copy phasefrac[cond_multi_index, phase_idx] * G'' to [cond_multi_index, np.index_exp[var_idx:last_var_idx, var_idx:last_var_idx]]
   Phasefrac/sitefrac derivatives
     Copy G' (clipped) to [cond_multi_index, last_var_idx+phase_idx, var_idx:last_var_idx] and [cond_multi_index, var_idx:last_var_idx, last_var_idx+phase_idx]
   Phasefrac/phasefrac derivative
      [cond_multi_index, last_var_idx+phase_idx, last_var_idx+phase_idx] = 0 (need to do this since we initialize as identity matrix)
    Constraint Jacobian:
    Need to track constraint_offset for each flattened condition multi_index (each participating phase increases it)
    Initialize constraint_jac as zeros of shape (*conds, max_phases, max_constraints, max_num_vars)
    Initialize l_constraints as zeros of shape(*conds, max_phases, max_constraints)
    Reset var_idx = first_var_idx
    Site fraction balances
      for idx in range(len(dbf.phases[name].sublattices)): active_in_subl = set(dbf.phases[name].constituents[idx]).intersection(comps)
      constraint_jac[cond_multi_index, phase_idx, constraint_offset, var_idx:last_var_idx+len(active_in_subl)] = 1
      l_constraints[cond_multi_index, phase_idx, constraint_offset] = -(np.sum(site_fracs[..., var_idx:var_idx + len(active_in_subl)], axis=-1) - 1)
      constraint_offset += 1
    Mass balances
      Basically the same as the existing code, with some indexing differences (remember max_phases dimension)
    Copy -constraint_jac.swapaxes(-2,-1).sum(axis=-3) to l_hessian[cond_multi_index, phase_idx, :max_num_vars, max_num_vars:]
    Copy constraint_jac.sum(axis=-3) to l_hessian[cond_multi_index, phase_idx, max_num_vars:, :max_num_vars]
    Gradient term:
    Reset var_idx = first_var_idx
    Initialize as b vector shape
    Copy -phasefrac[cond_multi_index, phase_idx] * G' (clipped) to [cond_multi_index, phase_idx, var_idx:var_idx+phase_dof[name]]
    Copy -G to [cond_multi_index, phase_idx, max_num_vars-max_phases+phase_idx]
    Copy np.einsum('...i, ...i', constraint_jac.swapaxes(-2,-1), l_multipliers[cond_multi_index, None, :]) to [cond_multi_index, phase_idx, :max_num_vars]
    Copy l_constraints to [cond_multi_index, phase_idx, max_num_vars:]

  Compute step = np.linalg.solve(l_hessian.sum(axis=-3), gradient_term.sum(axis=-2))
  step should have the shape (*conds, max_num_vars+num_constraints)
  Now, about the backtracking line search
  This is moderately annoying because we're going to need to compute GM for each condition set for different alpha's
  Initialize alpha as float ones of shape (*conds)

Blockers for 0.3 release

  • Rename equilibrium grid_opts kwarg to calc_opts
  • Add colored ZPF lines back to rewritten binplot (defer smooth interpolated lines to 0.4)
  • Add output kwarg support to equilibrium for lists of properties
  • Try to figure out cause of "singular matrix" error for one particular case in equilibrium (if cannot, defer to 0.3.x point release)
  • Regenerate API documentation because there have been a lot of breakages
  • Regenerate examples in documentation -- perhaps add one on degree-of-ordering calculation

Progress indicator for calculations

It would be nice to have some kind of indicator that a calculation was proceeding so that the user does not become concerned that it's stuck. It could be as easy as exposing some of the debugging output to the user, but I'm not sure how to write to stdout in an efficient way without interfering too much with the multi-threaded stuff happening in numpy, and also in a way that works across platforms and whether we're in the Notebook or in console mode.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.