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

Species support #161

Merged
merged 56 commits into from Mar 18, 2018
Merged

Species support #161

merged 56 commits into from Mar 18, 2018

Conversation

@richardotis
Copy link
Collaborator

@richardotis richardotis commented Feb 26, 2018

This PR adds support for charged species and associates in pycalphad. Several changes were made to the core. In particular,

  • A new Species object was added, which now replaces strings for the majority of species representation within the code (Database and Model in particular). Where we specifically refer to pure elements, we pull them out from the components list using a snippet (which should probably become a utility function).
  • Mass balance is now computed using a JIT-compiled function. This facilitated support for the new, more complex mass balance conditions for ionic liquids and associates. This also helps set the stage for having user-specified, JIT-compiled constraints (in a future PR).
  • To avoid further diverging code paths, CompiledModel is removed. We take a performance hit for this, but many of the original reasons this class was introduced are now moot due to JIT performance improvements. Also, CompiledModel performs poorly with complex phase models.

This PR turns PhaseRecord and Problem into a bit of a mess (dangling obsolete attributes and member functions, array allocation in the hot path, for example), but I wanted to make the minimal set of changes necessary to add the functionality since I intend to take another pass at these objects in the advanced constraints implementation.

Closes #110 and fixes #157

richardotis and others added 30 commits Jan 5, 2018
… breaking out species-specific mole fraction contributions
richardotis and others added 2 commits Mar 5, 2018
…ated code paths which are not taken. This makes wrapping constituent strings in Database.add_parameter() with Species() much cleaner.
Copy link
Collaborator

@bocklund bocklund left a comment

This finishes up the more formal 'code review'.

  • Hessian stuff is dead code with Ipopt, right? Since usage is scattered and it doesn't interfere with the current code, I'm ok with opening an issue to clean this up later.
  • What happens if users define conditions in terms of species?
  • What happens if we try to plot where a species is a dependent component in conditions?

Some things I still want to do:

  • run some test calculations in the BMG system
  • evaluate look and feel of returned Dataset objects
  • test some edge cases for different specified components/phases (e.g. phase has net charge with active components, such as (AL)(O-2) where there cannot be charge balance - TC suspends these phases by default)

I want to look for some edge cases and spend some time trying to break things, but other wise this looks good to me.

Loading

@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 5, 2018

  1. Yes, it's dead code.
  2. It might work?
  3. It should work.

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 6, 2018

On 2:

lih_tdb = """$ HLI
$
$ TDB-file for the thermodynamic assessment of the H-Li system
$
$--------------------------------------------------------------------------
$ 2008.12.11
$ 
$ TDB file for PANDAT created by M.Palumbo, T.Abe and K.Hashimoto
$
$ Thermodynamics Modeling Group , National Institute for 
$ Materials Science. 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan
$ e-mail: abe.taichi @nims.go.jp
$ Copyright (C) NIMS 2016
$
$
$ PARAMETERS ARE TAKEN FROM 
$   A. D. Pelton, Z. Metallkd. 84 (1993) p. 767-772
$
$ REMARKS: Assessment of the system for H content <= 50%at
$          A conversion in the liquid interaction parameters has been 
$          applied in order to use the standard Redlich-Kister polynomials
$          with respect to the polynomial form used in the original work
$
$ -------------------------------------------------------------------------
$
 ELEMENT /-   ELECTRON_GAS              0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT VA   VACUUM                    0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT H    1/2_MOLE_H2(G)            1.0079E+00  4.2340E+03  6.5285E+01!
 ELEMENT LI   BCC_A2                    6.9410E+00  4.6233E+03  2.9096E+01! 

 SPECIES H1LI1        H1LI1!
$
$ -------------------------------------------------------------------------
$
 FUNCTION GHSERHH  298.15   -4761.48705+39.26369395*T-15.678535*T*LN(T)
                            +.00137949625*T**2-3.731953335E-07*T**3
                            +28291.15*T**(-1);                 1000   Y
                            +90.054332-7.8064128*T-8.924285*T*LN(T)
                            -.00292084*T**2+1.573093335E-07*T**3
                             -640018*T**(-1);                  2100   Y 
                            -9420.08315+46.15601275*T-16.02541*T*LN(T)
                            -.00053641175*T**2+0.571408915E-08*T**3
                            +1780501.25*T**(-1);               6000   N !
 FUNCTION GHSERLI  200.00   -10583.817+217.636843*T-38.940488*T*LN(T)
                            +.035466931*T**2-1.9869816E-05*T**3
                            +159994*T**(-1);                   453.69 Y 
                            -565846.534+10668.0378*T-1722.32564*T*LN(T)
                            +2.28509159*T**2-5.77968153E-04*T**3
                            +34264368*T**(-1);                 500  Y 
                            -9062.949+179.277549*T-31.2283718*T*LN(T)
                            +.002633221*T**2-4.38058E-07*T**3
                            -102387*T**(-1);                   3000 N !
 FUNCTION GLILIQ   200.00   +2700.405-5.794896*T+GHSERLI#;     250  Y 
                            +577841.704-11029.6437*T+1783.83417*T*LN(T)
                            -2.46729546*T**2+6.41837109E-04*T**3
                            -34823613*T**(-1)+GHSERLI#;        453.69 Y 
                            +3005.582-6.624568*T+GHSERLI#;     3000   N !
 FUNCTION GLIFCC   200.00   -108+1.3*T+GHSERLI#;               3000   N !
 FUNCTION GLIHCP   200.00   -154+2*T+GHSERLI#;                 3000   N !
 FUNCTION GLIH     298.15   -112438+238.089*T-20.260*T*LN(T)
                            +.0035525*T**2-3.257E-07*T**3
                            +4.450E-05*T**(-1)+GLILIQ#+GHSERHH#;
                                                               3000   N !

 FUNCTION UN_ASS   298.15   0; 				        300   N !
 
 TYPE_DEFINITION % SEQ *!
 DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
 DEFAULT_COMMAND DEF_SYS_ELEMENT VA !
$
$ -------------------------------------------------------------------------
$
 PHASE LIQUID:L %  1  1.0  !
  CONSTITUENT LIQUID:L :LI,H1LI1,H :  !
  PARAMETER G(LIQUID,LI;0)        298.15  +GLILIQ#;    	   3000 N !
  PARAMETER G(LIQUID,H1LI1;0)     298.15  +GLIH#+25440-26.322*T; 3000 N !
  PARAMETER G(LIQUID,H1LI1,LI;0)  298.15  29049.75-8.7273*T;     3000 N !
  PARAMETER G(LIQUID,H1LI1,LI;1)  298.15  13803.5-6.1815*T;      3000 N !
  PARAMETER G(LIQUID,H1LI1,LI;2)  298.15  2640.75;               3000 N !

 PHASE BCC_A2  %  2 1   3 !
  CONSTITUENT BCC_A2  :LI% : VA% :  !
  PARAMETER G(BCC_A2,LI:VA;0)     298.15  +GHSERLI#;  	         3000 N !
$
$ Compound phases ----------------------------------------------------------
$
 PHASE LIH  %  2  1.0 1.0  !
  CONSTITUENT LIH  :LI :H :  !
  PARAMETER G(LIH,LI:H;0)         298.15  GLIH#;                 3000 N !
$
$------------------------------------------------------------ END OF LINE
 

"""

from pycalphad import equilibrium, variables as v, Database
dbf = Database(lih_tdb)
comps = ['H', 'H1LI1', 'VA', 'LI']
phases = list(dbf.phases.keys())
conds = {v.P:101325, v.X('LI'): (0, 1, 0.1), v.T: (500, 3000, 20)}
print(dbf.species)
eq = equilibrium(dbf, comps, phases, {v.P: 101325, v.T: 300, v.X('H1L1'): 0.5})

gives me

{Species('H', 'H1'), Species('/-', '/-1'), Species('H1LI1', 'H1.0LI1.0'), Species('VA', 'VA1'), Species('LI', 'LI1')}
---------------------------------------------------------------------------
ConditionError                            Traceback (most recent call last)
<ipython-input-49-eddc0e08b05e> in <module>()
      5 conds = {v.P:101325, v.X('LI'): (0, 1, 0.1), v.T: (500, 3000, 20)}
      6 print(dbf.species)
----> 7 eq = equilibrium(dbf, comps, phases, {v.P: 101325, v.T: 300, v.X('H1L1'): 0.5})

~/Projects/pycalphad/pycalphad/core/equilibrium.py in equilibrium(dbf, comps, phases, conditions, output, model, verbose, broadcast, calc_opts, scheduler, parameters, **kwargs)
    235     for cond in conds.keys():
    236         if isinstance(cond, (v.Composition, v.ChemicalPotential)) and cond.species not in comps:
--> 237             raise ConditionError('{} refers to non-existent component'.format(cond))
    238     str_conds = OrderedDict((str(key), value) for key, value in conds.items())
    239     num_calcs = np.prod([len(i) for i in str_conds.values()])

ConditionError: X_H1L1 refers to non-existent component

Loading

@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 6, 2018

I am on the fence about whether we should fix this now or wait until advanced conditions support.

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 6, 2018

Loading

@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 6, 2018

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 7, 2018

Compilation times with SymPy are not great for point calculations, but for repeated calculations or multiple condition mapping (phase diagrams) I think the compilation pays for itself quickly. As a demo tool, it's probably not ideal to switch between databases too much without being willing to talk through the compilation time.

I'm sure there are several chunk optimizations to be made, but even just parallelizing compiling across phases could be an easy win with a significant performance bump.

Loading

@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 8, 2018

Is parallel compile a blocker for this?

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 10, 2018

comps can not always be passed as string names to calculate. comps = ['ALO3/2'] will fail because possible_comps = {v.Species(x) for x in comps} cannot parse out a chemical formula and v.Species(x) returns None. Possible fix is to try to look up string compositions in a species dictionary from the database of {species_name: Species_obj}?

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 10, 2018

Are the kind of DataArrays you showed me OOB going to make it into this?

<xarray.DataArray 'Y' (P: 1, T: 1, points: 1, internal_dof: 5)>
array([[[[ 0.33333333,  0.33333333,  0.33333333,  1.        ,  0.        ]]]])
Coordinates:
  * T             (T) float64 300.0
  * P             (P) float64 1.013e+05
  * points        (points) int64 0
  * internal_dof  (internal_dof) int64 0 1 2 3 4

How should I inspect the site ratios of species-containing phases in calculate and equilibrium?

Loading

…to support species whose names do not parse to the chemical formula
@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 10, 2018

I'm thinking to defer those issues to advanced conditions, since it's where I'm also going to add support for things like weight fractions. Today I would resolve those issues (fraction of a species, site ratios of ionic liquids) by defining a custom property for the quantity you want.

Loading

bocklund added 2 commits Mar 10, 2018
Fixes error where eqplot tries to index Dataset components with Species objects
This is a performance optimization based on profiling where we can expect a modest speedup in tight equilibrium loops when callables are passed.
Since Hessian callables are not generated, hess_callable_dict.get(name, False) would always return None, triggering the if statement to always evaluate to True.
@richardotis
Copy link
Collaborator Author

@richardotis richardotis commented Mar 16, 2018

Is the species string parsing issue mentioned above fixed since 7562a86?

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 17, 2018

Yes, 7562a86 fixed parsing species names that are not chemical compositions and should close #157

Loading

@bocklund
Copy link
Collaborator

@bocklund bocklund commented Mar 18, 2018

This is ready to merge IMO

Loading

@richardotis richardotis merged commit 259184c into develop Mar 18, 2018
5 checks passed
Loading
@richardotis richardotis moved this from Review to Done in Roadmap Aug 18, 2018
@richardotis richardotis deleted the species-support branch Nov 14, 2018
bocklund added a commit to bocklund/pycalphad that referenced this issue Aug 17, 2021
…d#161)

Fixes a bug introduced by pycalphad#158 where GHSER data was set to zero. Includes a test.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
No open projects
Linked issues

Successfully merging this pull request may close these issues.

2 participants