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

Chemical potentials calculated by pycalphad are discontinuous #43

Closed
yangshenglan opened this issue May 18, 2016 · 20 comments
Closed

Chemical potentials calculated by pycalphad are discontinuous #43

yangshenglan opened this issue May 18, 2016 · 20 comments
Labels
bug
Milestone

Comments

@yangshenglan
Copy link

@yangshenglan yangshenglan commented May 18, 2016

By means of pycalphad, the curved surface of the chemical potentials calculated are discontinuous. As is shown in the figure1, this is the scatter diagram of the chemical potentials of component Ni of the Al-Ni-Cr ternary system at 1273K. In the figure1, some of the chemical potentials change dramatically in different compositions.
figure1
Selecting a composition corresponding to the dramatically variational chemical potential from figure1, calculate the chemical potentials of different compositions around that composition. As is shown figure2, this is the diagram of the chemical potentials of component Ni of compositions of the area.
figure2
How to solve the problem that the chemical potentials are discontinuous calculated by pycalphad?

Chemical potentials calculated by pycalphad are dis.docx

@yangshenglan yangshenglan changed the title Chemial potentials calculated by pycalphad are discontinuous Chemical potentials calculated by pycalphad are discontinuous May 18, 2016
@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 18, 2016

Can you verify that this problem exists in the latest release of pycalphad, 0.3.5?

@zhongjingjogy
Copy link

@zhongjingjogy zhongjingjogy commented May 18, 2016

Hi Otis,

It seems that we are working with the 0.3.5., and I installed this version for her the other day.

>>> pycalphad._version.get_versions()
{'version': '0.3.5-py2.7.egg', 'full': ''}

Best wishes!

@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 18, 2016

Can you provide the TDB and Python file used to perform this calculation so I can try to reproduce the bug?

As a workaround, you can try passing calc_opts={'pdens': 1000} to equilibrium()

@richardotis richardotis added the bug label May 18, 2016
@yangshenglan yangshenglan reopened this May 19, 2016
@yangshenglan
Copy link
Author

@yangshenglan yangshenglan commented May 19, 2016

TDB


$ Database file written 2016- 5-19
$ From database: TTNI5                   
 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 AL   FCC_A1                    2.6982E+01  4.5400E+03  2.8300E+01!
 ELEMENT CR   BCC_A2                    5.1996E+01  4.0500E+03  2.3543E+01!
 ELEMENT NI   FCC_A1                    5.8690E+01  4.7870E+03  2.9796E+01!


 FUNCTION GHSERAL   298.14 -7976.15+137.093038*T-24.3671976*T*LN(T)
     -.001884662*T**2-8.77664E-07*T**3+74092*T**(-1); 700 Y
      -11276.24+223.048446*T-38.5844296*T*LN(T)+.018531982*T**2
     -5.764227E-06*T**3+74092*T**(-1); 933.47 Y
      -11278.378+188.684153*T-31.748192*T*LN(T)-1.231E+28*T**(-9); 2900 N !
 FUNCTION ZERO      298.15 +0.0; 6000 N !
 FUNCTION GFCCCR    298.15 +7284+.163*T+GHSERCR#; 6000 N !
 FUNCTION GHSERNI   298.14 -5179.159+117.854*T-22.096*T*LN(T)-.0048407*T**2; 
     1728 Y
      -27840.655+279.135*T-43.1*T*LN(T)+1.12754E+31*T**(-9); 3000 N !
 FUNCTION GHSERCR   298.14 -8856.94+157.48*T-26.908*T*LN(T)+.00189435*T**2
     -1.47721E-06*T**3+139250*T**(-1); 2180 Y
      -34869.344+344.18*T-50*T*LN(T)-2.88526E+32*T**(-9); 6000 N !
 FUNCTION UN_ASS 298.15 +0; 300 N !

 TYPE_DEFINITION % SEQ *!
 DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
 DEFAULT_COMMAND DEF_SYS_ELEMENT VA /- !


 TYPE_DEFINITION & GES A_P_D FCC_A1 MAGNETIC  -3.0    2.80000E-01 !
 PHASE FCC_A1  %&  2 1   1 !
    CONSTITUENT FCC_A1  :AL,CR,NI% : VA% :  !

   PARAMETER G(FCC_A1,AL:VA;0)            298.15 +GHSERAL#; 6000 N REF0 !
   PARAMETER TC(FCC_A1,AL:VA;0)           298.15 +ZERO#; 6000 N REF0 !
   PARAMETER BMAGN(FCC_A1,AL:VA;0)        298.15 +ZERO#; 6000 N REF0 !
   PARAMETER G(FCC_A1,CR:VA;0)            298.15 +GFCCCR#; 6000 N REF0 !
   PARAMETER TC(FCC_A1,CR:VA;0)           298.15 -1109; 6000 N REF0 !
   PARAMETER BMAGN(FCC_A1,CR:VA;0)        298.15 -2.46; 6000 N REF0 !
   PARAMETER G(FCC_A1,NI:VA;0)            298.15 +GHSERNI#; 6000 N REF0 !
   PARAMETER TC(FCC_A1,NI:VA;0)           298.15 +633; 6000 N REF0 !
   PARAMETER BMAGN(FCC_A1,NI:VA;0)        298.15 +.52; 6000 N REF0 !
   PARAMETER G(FCC_A1,AL,CR:VA;0)         298.15 -45900+6*T; 6000 N REF0 !
   PARAMETER G(FCC_A1,AL,CR,NI:VA;0)      298.15 +30500-5*T; 6000 N REF0 !
   PARAMETER G(FCC_A1,AL,NI:VA;0)         298.15 -173950+19.35*T; 6000 N 
  REF0 !
   PARAMETER G(FCC_A1,AL,NI:VA;1)         298.15 +28500; 6000 N REF0 !
   PARAMETER G(FCC_A1,AL,NI:VA;2)         298.15 +47000; 6000 N REF0 !
   PARAMETER G(FCC_A1,CR,NI:VA;0)         298.15 +8030-12.8801*T; 6000 N 
  REF0 !
   PARAMETER G(FCC_A1,CR,NI:VA;1)         298.15 +33080-16.0362*T; 6000 N 
  REF0 !
   PARAMETER TC(FCC_A1,CR,NI:VA;0)        298.15 -3605; 6000 N REF0 !
   PARAMETER BMAGN(FCC_A1,CR,NI:VA;0)     298.15 -1.91; 6000 N REF0 !


 TYPE_DEFINITION ' GES A_P_D GAMMA_PRIME MAGNETIC  -3.0    2.80000E-01 !
 PHASE GAMMA_PRIME  %'  2 .75   .25 !
    CONSTITUENT GAMMA_PRIME  :AL,CR,NI% : AL%,CR,NI :  !

   PARAMETER G(GAMMA_PRIME,AL:AL;0)       298.15 +GHSERAL#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR:AL;0)       298.15 -8606+1.125*T+.75*GFCCCR#
  +.25*GHSERAL#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:AL;0)       298.15 -39860+3.175*T+.75*GHSERNI#
  +.25*GHSERAL#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL:CR;0)       298.15 -8606+1.125*T+.75*GHSERAL#
  +.25*GFCCCR#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR:CR;0)       298.15 +GFCCCR#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:CR;0)       298.15 -1915-.91*T+.75*GHSERNI#
  +.25*GFCCCR#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL:NI;0)       298.15 -35000+5*T+.75*GHSERAL#
  +.25*GHSERNI#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR:NI;0)       298.15 +10000-3.92*T+.75*GFCCCR#
  +.25*GHSERNI#; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:NI;0)       298.15 +GHSERNI#; 6000 N REF0 !
   PARAMETER TC(GAMMA_PRIME,NI:NI;0)      298.15 +633; 6000 N REF0 !
   PARAMETER BMAGN(GAMMA_PRIME,NI:NI;0)   298.15 +.52; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL,NI:AL;0)    298.15 -59500+20*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL,NI:AL;1)    298.15 +53350; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL:AL,NI;0)    298.15 +10000; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL:AL,CR;0)    298.15 -2869+.375*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR:AL,NI;0)    298.15 +10000; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR:AL,CR;0)    298.15 -2869+.375*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:AL,NI;0)    298.15 +9100-1.5*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:AL,NI;1)    298.15 -5400; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,NI:AL,CR;0)    298.15 -6250+T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL,NI:CR;0)    298.15 -60000+20*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL,NI:NI;0)    298.15 -50000; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,AL,CR:*;0)     298.15 -25819+3.375*T; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,*:CR,NI;0)     298.15 +3000; 6000 N REF0 !
   PARAMETER G(GAMMA_PRIME,CR,NI:*;0)     298.15 -1500; 6000 N REF0 !

 LIST_OF_REFERENCES
 NUMBER  SOURCE
  ! 

python file 1

import matplotlib.pyplot as plt
from pycalphad import equilibrium
from pycalphad import Database, Model
import pycalphad.variables as v
import numpy as np

db = Database('NiAlCr.tdb')
my_phases = ['FCC_A1', 'GAMMA_PRIME']

xal = np.linspace(0.12, 0.13, 50)
xcr = np.linspace(0.15, 0.1846, 100)
m = []
a = []
b = []
c = []
d = []
e = []
for i in range(50):
    for j in range(100):
        eq = equilibrium(db, ['AL', 'NI', 'CR', 'VA'], my_phases, {v.X('AL'): xal[i], v.X('CR'): xcr[j], v.T: 1273, v.P: 101325}, calc_opts={'pdens': 1000})
        m.append(eq)
        a.append(eq.get("X"))
        b.append(eq.get("Y"))
        c.append(eq.get("MU"))
        d.append(eq.get("GM"))
        e.append(eq.get("NP"))

f = open('equilibrium_data1.txt','w')
print >>f,m
print >>f,a
print >>f,b
print >>f,c
print >>f,d
print >>f,e
f.close()

python file 2

import matplotlib.pyplot as plt
from pycalphad import equilibrium
from pycalphad import Database, Model
import pycalphad.variables as v
import numpy as np

db = Database('NiAlCr.tdb')
my_phases = ['FCC_A1', 'GAMMA_PRIME']

xal = np.linspace(0.1211, 0.1233, 50)
xcr = np.linspace(0.15, 0.1788, 100)
m = []
a = []
b = []
c = []
d = []
e = []
for i in range(50):
    for j in range(100):
        eq = equilibrium(db, ['AL', 'NI', 'CR', 'VA'], my_phases, {v.X('AL'): xal[i], v.X('CR'): xcr[j], v.T: 1273, v.P: 101325}, calc_opts={'pdens': 1000})
        m.append(eq)
        a.append(eq.get("X"))
        b.append(eq.get("Y"))
        c.append(eq.get("MU"))
        d.append(eq.get("GM"))
        e.append(eq.get("NP"))
#
f = open('equilibrium_data2.txt','w')
print >>f,m
print >>f,a
print >>f,b
print >>f,c
print >>f,d
print >>f,e
f.close()
@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 19, 2016

Thanks. I will take a look at this problem this weekend.

@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 21, 2016

I can't reproduce this on Linux with Python 2.7 or 3.5, using either pycalphad 0.3.5 or the pycalphad develop branch. I cut your example down a bit to reduce the computation time. (The next version of pycalphad will be 100x faster for mapping.)

import matplotlib.pyplot as plt
from pycalphad import equilibrium
from pycalphad import Database, Model
import pycalphad.variables as v
import numpy as np

db = Database('NiAlCr.tdb')
my_phases = ['FCC_A1', 'GAMMA_PRIME']

xcr = np.linspace(0.17, 0.18, 10)

eq = equilibrium(db, ['AL', 'NI', 'CR', 'VA'], my_phases,
                 {v.X('AL'): .12, v.X('CR'): xcr, v.T: 1273, v.P: 101325})
print(eq)
fig = plt.figure(figsize=(12,12))
ax = fig.gca()
X = eq.X_CR.values
ax.set_xlabel('X(CR)')
Y = eq.MU.sel(component='AL')
ax.scatter(X, Y)
plt.show()

image

@yangshenglan
Copy link
Author

@yangshenglan yangshenglan commented May 27, 2016

The same to your example, I get the result that chemial potentials of component Al are continuous. However, chemical potentials of component Ni are discontinuous. Could you give the example about calculating chemical potentials of component Ni al the same composition.
xal = 0.12 xcr = np.linspace(0.17, 0.18, 10)

@zhongjingjogy
Copy link

@zhongjingjogy zhongjingjogy commented May 27, 2016

Hi,

I'm not sure whether we have make mistakes about the pycalphad version. So I ran some new tests anyway, as I have reconfirmed that a 0.3.5 version had been installed.
Codes are presented as following, which is the similiar to the previous one, but I varied the composition of the Al from 0.15 to 0.1846.

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

db = Database('NiAlCr.tdb')
my_phases = ['FCC_A1', 'GAMMA_PRIME']

x = np.empty((50, 100))
y = np.empty((50, 100))
mu_al, mu_cr, mu_ni = np.empty((50, 100)), np.empty((50, 100)), np.empty((50, 100))
gm = np.empty((50, 100))

xal = np.linspace(0.12, 0.13, 50)
xcr = np.linspace(0.15, 0.1846, 100)

for i in range(50):
    xcr = np.linspace(0.15, 0.1846, 100)
    eq = equilibrium(db, ['AL', 'NI', 'CR', 'VA'], my_phases,
                 {v.X('AL'): xal[i], v.X('CR'): xcr, v.T: 1273, v.P: 101325})
    print eq.X_AL
    print eq.MU.sel(component='AL')
    print eq.MU.sel(component='CR')
    print eq.MU.sel(component='NI')

ps: there are too much brakets in the output result, and I just manual came across the result manually rather than walking through the brakets with some iterations.

I came across the outputs the other day in order to figure out whether there are anything tricky. And I did observe some discontinous when the composition of Alumimum is 0.1241. Here follows the figure.
1241-all

Maybe the minimizer has been performed perfectly enough, and we are just stressing too much on the precision, which is not necessary.

What's more, there might a misibility gap in the Gibbs energy of FCC_A1 and I have no idea whether it bothers the calculation of the equilibirum or not.

Thanks!

@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 27, 2016

Thanks again for the details. I'll try to debug exactly this case and get back to you on the next 7-10 days.

@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 31, 2016

I can reproduce this error now. It appears to be related to the global minimization point density. If I increase pdens from default 100 to 200, the discontinuity goes away.

@richardotis
Copy link
Collaborator

@richardotis richardotis commented May 31, 2016

It appears to be a convergence failure in that composition range. I will look into this further.

@richardotis
Copy link
Collaborator

@richardotis richardotis commented Jun 8, 2016

@zhongjingjogy I want you to know I haven't forgotten about this and am still working on it. My PhD defense is coming up next week, so things have been very busy. I will probably not be able to release a fix for this until after next week, but I do intend to fix it.

@zhongjingjogy
Copy link

@zhongjingjogy zhongjingjogy commented Jun 16, 2016

That's very kind of you.

We might have figure out that the Gibbs energy of the gamma prime may have equivalent minima at different compositions(site fractions). It seems that one might have to take care of the casees that serveral minima exist in the Gibbs energy function of a single phase. So it might not be a problem of convergence, but a problem of chosing the wrong initials for the minimization algorithm.

Best wishes!

@richardotis
Copy link
Collaborator

@richardotis richardotis commented Jun 29, 2016

Quick update: I think I have this problem fixed in the develop branch. I should be able to release a new version in the next few days.

@richardotis
Copy link
Collaborator

@richardotis richardotis commented Jul 20, 2016

I've narrowed the problem down to an unexpected source: it looks like the code for removing composition sets is not working correctly. If I disable that code and turn up the global minimization point density, I get the correct answer. I'm working on a more robust solution now.

@zhongjingjogy
Copy link

@zhongjingjogy zhongjingjogy commented Jul 21, 2016

Thanks for your time! We have turn to the ThermoCalc for help at present, and we think pycalphad is an excellent open source project as long as it works in the right. I have been wondering why the TC do not provide any python interfaces at al, and I have to spend a half day on writing a script to parse the output from TC.

Best wishes!

@richardotis
Copy link
Collaborator

@richardotis richardotis commented Jul 21, 2016

Thanks @zhongjingjogy . It is not so easy to write a robust solver, so I have a lot of respect for the developers of Thermo-Calc. Once I've finally solved this issue I will let you know, and perhaps you can try pycalphad again for your next project.

richardotis added a commit that referenced this issue Aug 3, 2016
This still needs a way to add phases back to equilibrium if they were
missed in the initial starting point calculation. It is conceptually
pretty simple, but probably it will be easier if we wait after a code
cleanup.

In the mean time, we can increase the default point density as a
workaround.
@richardotis richardotis added this to the 0.4 milestone Aug 3, 2016
@richardotis
Copy link
Collaborator

@richardotis richardotis commented Aug 3, 2016

I have good news to report: I believe I have fixed this issue. Substantial portions of the solver core have been rewritten in the past week, and the result is that convergence is much more robust than before. I have also increased the default point density. Here are the fixed Al chemical potentials from your earlier test.
image

Thank you again for reporting. I have added your test case to pycalphad's test suite, so it should stay fixed permanently. There is still some work to do to add the ability for pycalphad to add new phases to the calculation during refinement, but with the solver core improvements and higher default point density, you should have no more problems calculating in this system. If you do have a problem, try increasing the value of pdens from 300 to something higher, like 1000. If that does not work, please file a new issue.

@zhongjingjogy I am closing this issue now since this is fixed in the develop branch, but I will send you a note when pycalphad 0.4 is released, so you can update.

@richardotis richardotis closed this Aug 3, 2016
richardotis added a commit that referenced this issue Aug 4, 2016
@richardotis
Copy link
Collaborator

@richardotis richardotis commented Aug 4, 2016

@zhongjingjogy pycalphad 0.4 has been released. Please update and let me know if it solves your problem.

@zhongjingjogy
Copy link

@zhongjingjogy zhongjingjogy commented Aug 4, 2016

Cheers! I am busy with some other business and I would like to have a try when I finish it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked pull requests

Successfully merging a pull request may close this issue.

None yet
3 participants
You can’t perform that action at this time.