In [2]:

    

from pycalphad.tests.datasets import *
from pycalphad.core.solver import InteriorPointSolver
from pycalphad import equilibrium, Database, Model, variables as v
import sympy

ALFE_DBF = Database(ALFE_TDB)
ISSUE43_DBF = Database('issue43.tdb')

class ProblemSaver(InteriorPointSolver):
    saved_problem = [None]
    def solve(self, prob):
        self.saved_problem[0] = prob
        self.verbose = True
        return super(ProblemSaver, self).solve(prob)

my_phases = ['FCC_A1', 'GAMMA_PRIME']
comps = ['AL', 'NI', 'CR', 'VA']
comps = sorted(comps)
conds = dict({v.T: 800, v.P: 101325, v.N: 1})



mod = Model(ISSUE43_DBF, comps, 'FCC_A1')

NP = sympy.Symbol('NP', real=True)
total_moles = sum([NP*mod.moles(c) for c in comps])
total_moles = NP

variables = [v.N, v.P, v.T] + mod.site_fractions + [NP]

mass_cons = [v.N, v.P, v.T]
mass_cons.extend(mod.get_internal_constraints())
mass_cons.extend(NP*mod.moles(c) for c in comps)

mass_jac = []
for cons in mass_cons:
    mass_jac.append([cons.diff(x) for x in variables])
    
energy_grad = [(total_moles*mod.GM).diff(x) for x in variables]


    

In [3]:

    

mass_cons


    

    
Out[3]:





[N,
 P,
 T,
 FCC_A10AL + FCC_A10CR + FCC_A10NI - 1,
 FCC_A11VA - 1,
 1.0*FCC_A10AL*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
 1.0*FCC_A10CR*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
 1.0*FCC_A10NI*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
 0]

In [4]:

    

mass_jac


    

    
Out[4]:





[[1, 0, 0, 0, 0, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0, 0],
 [0, 0, 1, 0, 0, 0, 0, 0],
 [0, 0, 0, 1, 1, 1, 0, 0],
 [0, 0, 0, 0, 0, 0, 1, 0],
 [0,
  0,
  0,
  -1.0*FCC_A10AL*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2 + 1.0*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
  -1.0*FCC_A10AL*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  -1.0*FCC_A10AL*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  0,
  1.0*FCC_A10AL/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)],
 [0,
  0,
  0,
  -1.0*FCC_A10CR*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  -1.0*FCC_A10CR*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2 + 1.0*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
  -1.0*FCC_A10CR*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  0,
  1.0*FCC_A10CR/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)],
 [0,
  0,
  0,
  -1.0*FCC_A10NI*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  -1.0*FCC_A10NI*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2,
  -1.0*FCC_A10NI*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)**2 + 1.0*NP/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI),
  0,
  1.0*FCC_A10NI/(1.0*FCC_A10AL + 1.0*FCC_A10CR + 1.0*FCC_A10NI)],
 [0, 0, 0, 0, 0, 0, 0, 0]]

In [5]:

    

A = sympy.Matrix(mass_jac).T.pinv()
x = A * sympy.Matrix(energy_grad)


    

    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-5-8724403478c9> in <module>
----> 1 A = sympy.Matrix(mass_jac).T.pinv()
      2 x = A * sympy.Matrix(energy_grad)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in pinv(self)
   4064                 return (AH * A).inv() * AH
   4065             else:
-> 4066                 return AH * (A * AH).inv()
   4067         except ValueError:
   4068             # Matrix is not full rank, so A*AH cannot be inverted.

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in inv(self, method, **kwargs)
   3192         if method is not None:
   3193             kwargs['method'] = method
-> 3194         return self._eval_inverse(**kwargs)
   3195 
   3196     def is_nilpotent(self):

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/dense.py in _eval_inverse(self, **kwargs)
    267         M = self.as_mutable()
    268         if method == "GE":
--> 269             rv = M.inverse_GE(iszerofunc=iszerofunc)
    270         elif method == "LU":
    271             rv = M.inverse_LU(iszerofunc=iszerofunc)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in inverse_GE(self, iszerofunc)
   3105 
   3106         big = Matrix.hstack(self.as_mutable(), Matrix.eye(self.rows))
-> 3107         red = big.rref(iszerofunc=iszerofunc, simplify=True)[0]
   3108         if any(iszerofunc(red[j, j]) for j in range(red.rows)):
   3109             raise ValueError("Matrix det == 0; not invertible.")

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in rref(self, iszerofunc, simplify, pivots, normalize_last)
    934         ret, pivot_cols = self._eval_rref(iszerofunc=iszerofunc,
    935                                           simpfunc=simpfunc,
--> 936                                           normalize_last=normalize_last)
    937         if pivots:
    938             ret = (ret, pivot_cols)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in _eval_rref(self, iszerofunc, simpfunc, normalize_last)
    589         reduced, pivot_cols, swaps = self._row_reduce(iszerofunc, simpfunc,
    590                                                       normalize_last, normalize=True,
--> 591                                                       zero_above=True)
    592         return reduced, pivot_cols
    593 

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in _row_reduce(self, iszerofunc, simpfunc, normalize_last, normalize, zero_above)
    705             pivot_offset, pivot_val, \
    706             assumed_nonzero, newly_determined = _find_reasonable_pivot(
--> 707                 get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
    708 
    709             # _find_reasonable_pivot may have simplified some things

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/matrices/matrices.py in _find_reasonable_pivot(col, iszerofunc, simpfunc)
   4757         if possible_zeros[i] is not None:
   4758             continue
-> 4759         simped = simpfunc(x)
   4760         is_zero = iszerofunc(simped)
   4761         if is_zero == True or is_zero == False:

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/simplify/simplify.py in simplify(expr, ratio, measure, rational, inverse)
    603 
    604     short = shorter(powsimp(expr, combine='exp', deep=True), powsimp(expr), expr)
--> 605     short = shorter(short, cancel(short))
    606     short = shorter(short, factor_terms(short), expand_power_exp(expand_mul(short)))
    607     if short.has(TrigonometricFunction, HyperbolicFunction, ExpBase):

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/polys/polytools.py in cancel(f, *gens, **args)
   6615 
   6616     try:
-> 6617         (F, G), opt = parallel_poly_from_expr((p, q), *gens, **args)
   6618     except PolificationFailed:
   6619         if not isinstance(f, (tuple, Tuple)):

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/polys/polytools.py in parallel_poly_from_expr(exprs, *gens, **args)
   4307     """Construct polynomials from expressions. """
   4308     opt = options.build_options(gens, args)
-> 4309     return _parallel_poly_from_expr(exprs, opt)
   4310 
   4311 

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/polys/polytools.py in _parallel_poly_from_expr(exprs, opt)
   4346 
   4347                 if opt.expand:
-> 4348                     expr = expr.expand()
   4349         else:
   4350             failed = True

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/cache.py in wrapper(*args, **kwargs)
     92             def wrapper(*args, **kwargs):
     93                 try:
---> 94                     retval = cfunc(*args, **kwargs)
     95                 except TypeError:
     96                     retval = func(*args, **kwargs)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/expr.py in expand(self, deep, modulus, power_base, power_exp, mul, log, multinomial, basic, **hints)
   3198             if hints.get('mul', False):
   3199                 expr, _ = Expr._expand_hint(
-> 3200                     expr, '_eval_expand_mul', deep=deep, **hints)
   3201             if hints.get('log', False):
   3202                 expr, _ = Expr._expand_hint(

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3122             sargs = []
   3123             for arg in expr.args:
-> 3124                 arg, arghit = Expr._expand_hint(arg, hint, **hints)
   3125                 hit |= arghit
   3126                 sargs.append(arg)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3130 
   3131         if hasattr(expr, hint):
-> 3132             newexpr = getattr(expr, hint)(**hints)
   3133             if newexpr != expr:
   3134                 return (newexpr, True)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/mul.py in _eval_expand_mul(self, **hints)
    877             if sums:
    878                 deep = hints.get("deep", False)
--> 879                 terms = self.func._expandsums(sums)
    880                 args = []
    881                 for term in terms:

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/mul.py in _expandsums(sums)
    842         right = Mul._expandsums(sums[L//2:])
    843 
--> 844         terms = [Mul(a, b) for a in left for b in right]
    845         added = Add(*terms)
    846         return Add.make_args(added)  # it may have collapsed down to one term

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/mul.py in <listcomp>(.0)
    842         right = Mul._expandsums(sums[L//2:])
    843 
--> 844         terms = [Mul(a, b) for a in left for b in right]
    845         added = Add(*terms)
    846         return Add.make_args(added)  # it may have collapsed down to one term

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/cache.py in wrapper(*args, **kwargs)
     92             def wrapper(*args, **kwargs):
     93                 try:
---> 94                     retval = cfunc(*args, **kwargs)
     95                 except TypeError:
     96                     retval = func(*args, **kwargs)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/operations.py in __new__(cls, *args, **options)
     45             return args[0]
     46 
---> 47         c_part, nc_part, order_symbols = cls.flatten(args)
     48         is_commutative = not nc_part
     49         obj = cls._from_args(c_part + nc_part, is_commutative)

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/mul.py in flatten(cls, seq)
    606 
    607         # order commutative part canonically
--> 608         _mulsort(c_part)
    609 
    610         # current code expects coeff to be always in slot-0

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/mul.py in _mulsort(args)
     32 def _mulsort(args):
     33     # in-place sorting of args
---> 34     args.sort(key=_args_sortkey)
     35 
     36 

~/anaconda3/envs/symengine/lib/python3.7/site-packages/sympy-1.4-py3.7.egg/sympy/core/basic.py in compare(self, other)
    206             l = Basic(*l) if isinstance(l, frozenset) else l
    207             r = Basic(*r) if isinstance(r, frozenset) else r
--> 208             if isinstance(l, Basic):
    209                 c = l.compare(r)
    210             else:

KeyboardInterrupt: 

In [5]:

    

from pycalphad.codegen.sympydiff_utils import build_functions

mu_a = build_functions(x[4], variables, include_grad=True, include_hess=False)
mu_b = build_functions(x[5], variables, include_grad=True, include_hess=False)
mu_c = build_functions(x[6], variables, include_grad=True, include_hess=False)

energy = build_functions(mod.GM, variables, include_grad=True)


    

In [6]:

    

print(x[4].free_symbols)


    

    




{TEST0A, TEST0B, NP, T, TEST0C}

In [7]:

    

mu_a.func([1, 1e5, 1000, 0.4, 0.6, 1e-12, 1e-6])


    

    
Out[7]:





array(-7618.14919886)

In [8]:

    

np.array(mu_a.grad([1, 1e5, 1000, 0.4, 0.6, 1e-12, 1])) - np.array(mu_b.grad([1, 1e5, 1000, 0.4, 0.6, 1e-12, 1]))


    

    
Out[8]:





array([ 0.00000000e+00,  0.00000000e+00, -3.37123964e+00,  2.07862500e+04,
       -1.38575000e+04,  9.23094935e-01, -9.44368765e-10])

In [9]:

    

from pycalphad.core.solver import InteriorPointSolver

class ProblemSaver(InteriorPointSolver):
    saved_problem = [None]
    def solve(self, prob):
        self.saved_problem[0] = prob
        self.verbose = True
        return super(ProblemSaver, self).solve(prob)

eq = equilibrium(dbf, ['A', 'B', 'C'], ['TEST'],
                     {v.MU('B'): -1000, v.X('A'): 0.1, v.T: 800, v.P: 101325}, solver=ProblemSaver())


    

    




Chemical Potentials [-15315.87500456  -1000.         -21480.14360388]
[0. 0. 0. 0. 0. 0. 0.]
[1.00000000e+00 1.01325000e+05 8.00000000e+02 1.00000000e-01
 8.60415585e-01 3.95844148e-02 1.00000000e+00]
Status: 0 b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'

In [10]:

    

ProblemSaver.saved_problem[0].jacobian([1, 1e5, 800, 0.1, 8.60415585e-01, 3.95844148e-2, 1.0])[-1]


    

    
Out[10]:





array([ 0.00000000e+00,  0.00000000e+00, -1.25000000e+00,  2.59387889e-09,
        7.73068284e+03,  1.93540473e-09, -8.77662387e-10])

In [11]:

    

mu_b.grad([1, 1e5, 800, 0.1, 8.60415585e-01, 3.95844148e-2, 1.0])


    

    
Out[11]:





[array(0.),
 array(0.),
 array(-1.25),
 array(-3.52429197e-12),
 array(7730.68284206),
 array(2.79669621e-11),
 array(2.33626452e-11)]

In [ ]: