```python
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]
```


```python
mass_cons
```




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




```python
mass_jac
```




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




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


```python
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)
```


```python
print(x[4].free_symbols)
```

    {TEST0A, TEST0B, NP, T, TEST0C}



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




    array(-7618.14919886)




```python
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]))
```




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




```python
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).'



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




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




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




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




```python

```

5cf77c30ee29c596e20106c58990818e7c1c955c
ipynb
Jupyter Notebook
SymbolicSolve.ipynb
richardotis/pycalphad-sandbox
43d8786eee8f279266497e9c5f4630d19c893092
MIT
2017-03-08T18:21:30.000Z
2017-03-08T18:21:30.000Z
SymbolicSolve.ipynb
richardotis/pycalphad-sandbox
43d8786eee8f279266497e9c5f4630d19c893092
MIT
SymbolicSolve.ipynb
richardotis/pycalphad-sandbox
43d8786eee8f279266497e9c5f4630d19c893092
MIT
2018-11-03T01:31:57.000Z
2018-11-03T01:31:57.000Z
Qwen/Qwen-72B
1. YES
2. YES