In [1]:

    

from pycalphad import equilibrium, Database, Model, variables as v
import sympy
import numpy as np

TDB = """
 ELEMENT A    GRAPHITE                   12.011     1054.0      5.7423 !
 ELEMENT B   BCC_A2                     55.847     4489.0     27.2797 !
 ELEMENT C   BCC_A2                     55.847     4489.0     27.2797 !
 TYPE_DEFINITION % SEQ * !
 PHASE TEST % 1 1 !
 CONSTITUENT TEST : A,B,C: !
"""
my_phases = ['TEST']
comps = ['A', 'B','C']
comps = sorted(comps)
conds = dict({v.T: 1000, v.P: 101325, v.N: 1})

dbf = Database(TDB)


mod = Model(dbf, ['A', 'B', 'C'], 'TEST')

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 [2]:

    

mass_cons


    

    
Out[2]:





[N,
 P,
 T,
 TEST0A + TEST0B + TEST0C - 1,
 1.0*TEST0A*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C),
 1.0*TEST0B*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C),
 1.0*TEST0C*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)]

In [3]:

    

mass_jac


    

    
Out[3]:





[[1, 0, 0, 0, 0, 0, 0],
 [0, 1, 0, 0, 0, 0, 0],
 [0, 0, 1, 0, 0, 0, 0],
 [0, 0, 0, 1, 1, 1, 0],
 [0,
  0,
  0,
  -1.0*TEST0A*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2 + 1.0*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C),
  -1.0*TEST0A*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  -1.0*TEST0A*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  1.0*TEST0A/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)],
 [0,
  0,
  0,
  -1.0*TEST0B*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  -1.0*TEST0B*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2 + 1.0*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C),
  -1.0*TEST0B*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  1.0*TEST0B/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)],
 [0,
  0,
  0,
  -1.0*TEST0C*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  -1.0*TEST0C*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2,
  -1.0*TEST0C*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)**2 + 1.0*NP/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C),
  1.0*TEST0C/(1.0*TEST0A + 1.0*TEST0B + 1.0*TEST0C)]]

In [4]:

    

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


    

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 [ ]: