Equilibrium fails to calculate, resulting in an IndexError

[bocklund]

The following code example fails to calculate with the given database. The database is known to work as tested in other software.

dbf = Database('/path/to/db.tdb')
phases = ['FCC'_A1, 'FCC_L12' ]
components = ["AL", "CO", "CR", "W", "VA"]
conditions = {"T": 1248, "P":101325, v.X("AL"):0.081, v.X("CR"): 0.020, v.X("W"):0.094}
equilibrium(dbf, components, phases, conditions)

The traceback for this is as follows

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-3-2193fcd054eb> in <module>()
     15 components = ["AL", "CO", "CR", "W", "VA"]
     16 conditions = {"T": 1248, "P":101325, v.X("AL"):0.081, v.X("CR"): 0.020, v.X("W"):0.094}
---> 17 equilibrium(dbf, components, phases, conditions)
     18 #equilibrium(dbf, ['AL', 'CO', 'CR', 'W'], phases, {v.X('CR'):0.05 v.X('CO'):0.77, v.X('W'):0.09, v.X('AL'):0.09, v.T:1250, v.P:101325})
     19 #binplot(db, ['AL', 'CR', 'VA'] , phases, {v.X('CR'):(0,1,0.02),

/Users/brandon/Projects/pycalphad/pycalphad/core/equilibrium.py in equilibrium(dbf, comps, phases, conditions, output, model, verbose, broadcast, calc_opts, scheduler, parameters, **kwargs)
    356     delayed(properties.attrs.__setitem__, pure=False)('created', datetime.utcnow())
    357     if scheduler is not None:
--> 358         properties = dask.compute(properties, get=scheduler)[0]
    359     return properties

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/base.py in compute(*args, **kwargs)
    176         dsk = merge(var.dask for var in variables)
    177     keys = [var._keys() for var in variables]
--> 178     results = get(dsk, keys, **kwargs)
    179 
    180     results_iter = iter(results)

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in get_sync(dsk, keys, **kwargs)
    548     kwargs.pop('num_workers', None)    # if num_workers present, remove it
    549     return get_async(apply_sync, 1, dsk, keys,
--> 550                      raise_on_exception=True, **kwargs)
    551 
    552 

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in get_async(apply_async, num_workers, dsk, result, cache, get_id, raise_on_exception, rerun_exceptions_locally, callbacks, dumps, loads, **kwargs)
    507 
    508             while state['ready'] and len(state['running']) < num_workers:
--> 509                 fire_task()
    510 
    511     except KeyboardInterrupt:

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in fire_task()
    474                         args=(key, dumps((dsk[key], data)),
    475                               dumps, loads, get_id, raise_on_exception),
--> 476                         callback=queue.put)
    477 
    478         # Seed initial tasks into the thread pool

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in apply_sync(func, args, kwds, callback)
    536 def apply_sync(func, args=(), kwds={}, callback=None):
    537     """ A naive synchronous version of apply_async """
--> 538     res = func(*args, **kwds)
    539     if callback is not None:
    540         callback(res)

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in execute_task(key, task_info, dumps, loads, get_id, raise_on_exception)
    266     try:
    267         task, data = loads(task_info)
--> 268         result = _execute_task(task, data)
    269         id = get_id()
    270         result = dumps((result, None, id))

/Users/brandon/anaconda3/envs/espei/lib/python3.5/site-packages/dask/async.py in _execute_task(arg, cache, dsk)
    247         func, args = arg[0], arg[1:]
    248         args2 = [_execute_task(a, cache) for a in args]
--> 249         return func(*args2)
    250     elif not ishashable(arg):
    251         return arg

/Users/brandon/Projects/pycalphad/pycalphad/core/eqsolver.pyx in pycalphad.core.eqsolver._solve_eq_at_conditions (pycalphad/core/eqsolver.c:17825)()

IndexError: index 4 is out of bounds for axis 6 with size 4

This was also confirmed using each phase individually and in another database and the same conditions, so it does not seem to be a problem with how any one specific model is represented.

This is on the develop branch, commit a518e1a

[richardotis]

I received the TDB out-of-band since the OP cannot distribute it publicly.

[richardotis]

With verbose=True

Calculation Backend: Compiled (autowrap)
Components: AL CO CR VA W
Phases: FCC_A1 [done]
('NEW_L_MULTIPLIERS', array([  7.72513013e+01,  -1.61106289e+04,   2.68604203e+03,
        -1.22397210e+04,  -1.44610437e+02,  -3.19581191e+04,
         1.92819422e+03,  -5.16770763e+03,  -1.49442772e+05,
        -6.46427102e+04,  -8.53929875e+04,  -6.58803572e+04]))
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-4-135e93af6fa7> in <module>()
      5 components = ["AL", "CO", "CR", "W", "VA"]
      6 conditions = {v.T: 1248, v.P:101325, v.X("AL"):0.081, v.X("CR"): 0.020, v.X("W"):0.094}
----> 7 equilibrium(dbf, components, phases, conditions, verbose=True)

/home/rotis/git/pycalphad/pycalphad/core/equilibrium.py in equilibrium(dbf, comps, phases, conditions, output, model, verbose, broadcast, calc_opts, scheduler, parameters, **kwargs)
    333         # Single-process job; don't create child processes
    334         properties = delayed(_solve_eq_at_conditions, pure=False)(dbf, comps, properties, phase_records,
--> 335                                                                  list(str_conds.keys()), verbose, diagnostic)
    336 
    337     # Compute equilibrium values of any additional user-specified properties

/home/rotis/git/pycalphad/pycalphad/core/eqsolver.pyx in pycalphad.core.eqsolver._solve_eq_at_conditions (pycalphad/core/eqsolver.c:18432)()

IndexError: index 4 is out of bounds for axis 6 with size 4

[richardotis]

It looks like the problem may be with some code that implicitly assumes that VA will always be the last component to appear alphabetically. The component index is incorrectly incremented. Working on a patch.

[richardotis]

Minimal test case which fails with the same error on develop:

def test_eq_issue62_last_component_not_va():
    """
    VA is not last when components are sorted alphabetically.
    """
    test_tdb = """
    ELEMENT VA   VACUUM                    0.0000E+00  0.0000E+00  0.0000E+00!
    ELEMENT AL   FCC_A1                    2.6982E+01  4.5773E+03  2.8322E+01!
    ELEMENT CO   HCP_A3                    5.8933E+01  4.7656E+03  3.0040E+00!
    ELEMENT CR   BCC_A2                    5.1996E+01  4.0500E+03  2.3560E+01!
    ELEMENT W    BCC_A2                    1.8385E+02  4.9700E+03  3.2620E+01!
    PHASE FCC_A1  %  2 1   1 !
    CONSTITUENT FCC_A1  :AL,CO,CR,W : VA% :  !
    """
    equilibrium(Database(test_tdb), ['AL', 'CO', 'CR', 'W', 'VA'], ['FCC_A1'],
                {"T": 1248, "P": 101325, v.X("AL"): 0.081, v.X("CR"): 0.020, v.X("W"): 0.094})