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.

Sign up for GitHub

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

Jump to bottom

still problem with solver #102

Closed
broshe opened this issue Jun 18, 2017 · 16 comments
Closed

still problem with solver #102

broshe opened this issue Jun 18, 2017 · 16 comments

Comments

@broshe
Copy link

broshe commented Jun 18, 2017

Hello Richard and Brandon,
I still find problems with the equilibrium solver.
See below some apparently simple calculations that give wrong results.
The calculations are done on the Al-Fe system on the Al-rich side.
At several temperatures, the calculation results in a wrong composition for the FCC_A1 phase.
Here is the script (after the script I give the results):
#!/usr/bin/env python3

-- coding: utf-8 --

"""
Created on Sun Jun 18 13:45:49 2017

@author: ebrosh
"""

from pycalphad import Database, equilibrium
import pycalphad.variables as v

db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')

print('conditions')
conditions={v.X('FE'): .0028156, v.T: 931.385693359375, v.P: 1.0e5}
print(conditions)
data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium')
print(data['Phase'].values.squeeze())
print('compositions of equilibrium phases')
print(data['X'].values.squeeze())

print('conditions')
conditions={v.X('FE'): .00460485, v.T: 929.933, v.P: 1.0e5}
print(conditions)
data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium')
print(data['Phase'].values.squeeze())
print('compositions of equilibrium phases')
print(data['X'].values.squeeze())

print('conditions')
conditions={v.X('FE'): .0047534, v.T: 929.933, v.P: 1.0e5}
print(conditions)
data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium')
print(data['Phase'].values.squeeze())
print('compositions of equilibrium phases')
print(data['X'].values.squeeze())

print('conditions')
conditions={v.X('FE'): .005477, v.T: 929.518, v.P: 1.0e5}
print(conditions)
data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium')
print(data['Phase'].values.squeeze())
print('compositions of equilibrium phases')
print(data['X'].values.squeeze())

this last equilibrium is correct

print('conditions')
conditions={v.X('FE'): .005477, v.T: 929, v.P: 1.0e5}
print(conditions)
data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium')
print(data['Phase'].values.squeeze())
print('compositions of equilibrium phases')
print(data['X'].values.squeeze())
print('note that the last one is correct')

Results of running the script:
conditions
{P: 100000.0, T: 931.385693359375, X_FE: 0.0028156}
phases in equilibrium
['FCC_A1' 'LIQUID']
compositions of equilibrium phases
[[ 1.00000000e+00 2.66513395e-11]
[ 9.97001429e-01 2.99857143e-03]]
conditions
{P: 100000.0, T: 929.933, X_FE: 0.00460485}
phases in equilibrium
['FCC_A1' 'LIQUID']
compositions of equilibrium phases
[[ 1.00000000e+00 2.75875971e-11]
[ 9.95001245e-01 4.99875457e-03]]
conditions
{P: 100000.0, T: 929.933, X_FE: 0.0047534}
phases in equilibrium
['FCC_A1' 'LIQUID']
compositions of equilibrium phases
[[ 1.00000000e+00 2.74369001e-11]
[ 9.95001245e-01 4.99875456e-03]]
conditions
{P: 100000.0, T: 929.518, X_FE: 0.005477}
phases in equilibrium
['FCC_A1' 'LIQUID']
compositions of equilibrium phases
[[ 1.00000000e+00 2.65337850e-11]
[ 9.94442085e-01 5.55791522e-03]]
conditions
{P: 100000.0, T: 929, X_FE: 0.005477}
phases in equilibrium
['FCC_A1' 'LIQUID']
compositions of equilibrium phases
[[ 9.99845081e-01 1.54919219e-04]
[ 9.93613993e-01 6.38600654e-03]]
note that the last one is correct

Can this be fixed ?

Best regards,
Eli
@bocklund
Copy link
Collaborator

Sorry, Eli, I'm a little confused on what exactly you are saying is incorrect. Are you saying all of those results are incorrect except for the last one?

What were you expecting to see?

Note that you can format your code nicely by wrapping it in three backticks, e.g.

```
eq = equilibrium(dbf, comps, phases, conditions)
```

would format to

eq = equilibrium(dbf, comps, phases, conditions)

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

Which version of pycalphad are you running?

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

bocklund commented Jun 18, 2017 •
edited

I am running on the HEAD of the develop branch (not much further than the 0.5.1 release) and on the Al-Fe_sei database straight from NIMS I have the following results that seem to be more in the ballpark (default pdens of 500):

from pycalphad import Database, equilibrium, variables as v
db_file = '/Users/brandon/Box Sync/databases/nims-databases/Al-Fe_sei.tdb'
dbf = Database(db_file)
comps = ['AL', 'FE', 'VA']
phases = ['FCC_A1', 'LIQUID']
conditions = {v.X('FE'): .0028156, v.T: 931.385693359375, v.P: 1.0e5}
eq = equilibrium(dbf, comps, phases, conditions)
eq.X
<xarray.DataArray 'X' (P: 1, T: 1, X_FE: 1, vertex: 2, component: 2)>
array([[[[[  9.99929563e-01,   7.04372289e-05],
          [  9.96935045e-01,   3.06495537e-03]]]]])
Coordinates:
  * X_FE       (X_FE) float64 0.002816
  * component  (component) <U2 'AL' 'FE'
  * T          (T) float64 931.4
  * P          (P) float64 1e+05
  * vertex     (vertex) int64 0 1

Then same conditions with pdens=2000 gives me the same answer that you were getting.

eq = equilibrium(dbf, comps, phases, conditions, calc_opts={'pdens': 2000})
<xarray.DataArray 'X' (P: 1, T: 1, X_FE: 1, vertex: 2, component: 2)>
array([[[[[  1.00000000e+00,   2.66513395e-11],
          [  9.97001429e-01,   2.99857143e-03]]]]])
Coordinates:
  * X_FE       (X_FE) float64 0.002816
  * component  (component) <U2 'AL' 'FE'
  * T          (T) float64 931.4
  * P          (P) float64 1e+05
  * vertex     (vertex) int64 0 1

Increasing the point density to 10000 goes back to the answer calculated by pdens=500:

eq = equilibrium(dbf, comps, phases, conditions, calc_opts={'pdens': 10000})
<xarray.DataArray 'X' (P: 1, T: 1, X_FE: 1, vertex: 2, component: 2)>
array([[[[[  9.99929563e-01,   7.04372289e-05],
          [  9.96935045e-01,   3.06495537e-03]]]]])
Coordinates:
  * X_FE       (X_FE) float64 0.002816
  * component  (component) <U2 'AL' 'FE'
  * T          (T) float64 931.4
  * P          (P) float64 1e+05
  * vertex     (vertex) int64 0 1

So it looks like that particular point density is giving some problems. Some coming changes by @richardotis should improve the solver further.

@bocklund
Copy link
Collaborator

Further, it looks like the 2000 point density isn't taking steps (verbose output):

Calculation Backend: Compiled (autowrap)
Components: AL FE VA
Phases: FCC_A1 LIQUID [done]
('NEW_L_MULTIPLIERS', array([    472.5425895 ,   -1829.33937232,    7065.81999672,
        -37723.13634075, -147658.20014453]))
('old_driving_force', -38032.669484192571)
(1.0, -38032.669495361733, 2.6818302023556839e-09)
('alpha', 1.0)
('Phases', [CompositionSet(FCC_A1, [  1.00000000e+00   2.66513395e-11]), CompositionSet(LIQUID, [ 0.99700143  0.00299857])])
('step', array([ -2.56502130e-11,   2.56513395e-11,  -1.44832580e-16,
         2.92931975e-06,  -2.92931976e-06,  -9.15504689e-04,
         9.15504689e-04]))
('Site fractions', array([  1.00000000e+00,   2.66513395e-11,   1.00000000e+00,
         9.97001429e-01,   2.99857143e-03]))
('Phase fractions', array([ 0.06102043,  0.93897957]))
('Chemical potentials', array([ -37723.13634075, -147658.20014453]))
('Chem pot progress', array([ -3.58490070e-06,  -6.61725708e-03]))
('Energy progress', 0.00028365789330564439)
('Driving force', 1.1037191143259406e-05)
No progress

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

bocklund commented Jun 18, 2017 •
edited

You could try a higher pdens that is not 2000. E.g. 4000 worked for me for the conditions

{T: 929.518, X_FE: 0.005477, P: 100000.0}

<xarray.DataArray 'X' (P: 1, T: 1, X_FE: 1, vertex: 2, component: 2)>
array([[[[[  9.998638e-01,   1.362262e-04],
          [  9.943198e-01,   5.680247e-03]]]]])
Coordinates:
  * P          (P) float64 1e+05
  * T          (T) float64 929.5
  * X_FE       (X_FE) float64 0.005477
  * vertex     (vertex) int64 0 1
  * component  (component) <U2 'AL' 'FE'

Indeed, increasing the point density should only be positive and it should never make your answer worse. My understanding is (@richardotis could explain better) that in the last example, with both a point density of 500 and 2000, your initial sampling of the energy surface is 'lucky' (unlucky in this case) in that the points that are sampled are the global minimum under the equilibrium conditions in the solver.

In both cases with the conditions {T: 929.518, X_FE: 0.005477, P: 100000.0} and pdens=500 or pdens=2000 you are getting convergence on the initial step.

I think this has to do with some trickiness of small solubilities, the slope of the free energy surface in that area, and the sampling resolution in that small region.

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

bocklund commented Jun 18, 2017 •
edited

The conditions are controlled in eqsolver.pyx. The conditions themselves are at the top of the file lines 14-21 and they are checked in lines 485-488.

There is currently no way to set these programmatically from equilibrium() or calculate() because they are hardcoded there. You are, of course, free to change them at your own risk.

In this case you can see that the driving force criteria is that any driving force below 1e-4 is converged and the points we have on the initial step from my comment above with pdens=2000 happens to give driving force of 1.103e-5 so it is 'converged' even though the actual answer is incorrect. The energy progress is also below the threshold.

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

.pyx files are compiled by Cython at build time. The files that end in *.so are the compiled versions of eqsolver.

In order to change the values you'll have to install pycalphad from source see the installation guide for installing development versions.

Every time you make a change in eqsolver.pyx, you'll have to recompile. The easiest way is to go to pycalphad's root directory (with the setup.py) and run python setup.py build_ext

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

@bocklund
Copy link
Collaborator

We will consider that suggestion, thank you!

Let us know if you need anything else. Ok to close this?

@broshe
Copy link
Author

broshe commented Jun 18, 2017 via email

bocklund added a commit to bocklund/pycalphad that referenced this issue Aug 17, 2021
@bocklund
MAINT: Remove automatically adding ideal exclusions (pycalphad#168)
f876b49 
The `add_ideal_exclusions` function was added by pycalphad#102 in June 2019 and was also deprecated at that time for removal in ESPEI 0.8. This function and its usages are now removed in favor of manually adding the exclusions or using tags.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@bocklund @broshe