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

Bug with custom model calculating entropy #452

Open
ianhbell opened this issue Jan 5, 2023 · 7 comments · May be fixed by #453
Open

Bug with custom model calculating entropy #452

ianhbell opened this issue Jan 5, 2023 · 7 comments · May be fixed by #453

Comments

@ianhbell
Copy link

ianhbell commented Jan 5, 2023

I am a newbie to pycalphad, but I believe I have run into a problem. I learned later on that calculation of the entropy was already available in the library, but in so doing, found a bug. If I calculate entropy from its definition based upon the Gibbs energy as $s=-\partial g/\partial T$, I always get 0 as the output, but no error message is returned. I derived my example from the viscosity calculation example which was very informative.

Code:

from tinydb import where
from pycalphad import Model, variables as v
from pycalphad import Database
from pycalphad import calculate
import pandas
import numpy as np

with open('unary50edt.tdb') as fp:
    dbf = Database(fp.read().upper())

class MyModel(Model):
    def build_phase(self, dbe):
        super(MyModel, self).build_phase(dbe)
        self.entropyval = self.build_entropy(dbe)

    def build_entropy(self, _):
        Gibbs_energy = self.GM
        entropy = -Gibbs_energy.diff(v.T)
        return entropy

for atom in ['Ni']:
    mod = MyModel(dbf, [atom], 'LIQUID')

    # NOTICE: we need to tell pycalphad about our model for this phase
    models = {'LIQUID': mod}

    o = []
    for temp in np.arange(950, 2000, 200):
        res = calculate(dbf, [atom], 'LIQUID', P=101325, T=temp, model=models, output='SM')
        sA = res.SM.values.squeeze()
        res = calculate(dbf, [atom], 'LIQUID', P=101325, T=temp, model=models, output='entropyval')
        sB = res.entropyval.values.squeeze()
        print(res)
        print(temp, sA, sB)
        
    pandas.DataFrame(o).to_csv(f'{atom}_CALPHAD_test.csv', index=False)

yielding

950 74.3563867319028 -0.0
1150 80.55645325671628 -0.0
1350 86.13576559192114 -0.0
1550 91.33373041483597 -0.0
1750 96.34571598414678 -0.0
1950 101.00972148214086 -0.0
@richardotis
Copy link
Collaborator

richardotis commented Jan 5, 2023

Good find! I suspect the root cause of the issue is this (taken from model.py, starting at line 255):

self.models = OrderedDict()
self.build_phase(dbe)

for name, value in self.models.items():
    # XXX: xreplace hack because SymEngine seems to let Symbols slip in somehow
    self.models[name] = self.symbol_replace(value, symbols).xreplace(v.supported_variables_in_databases)

Your custom build_phase is getting called after self.models is initialized, but before symbol replacement happens. Your call to self.GM uses self.models under the hood, so effectively you're accessing it before it's fully initialized. When you call .diff(v.T), I suspect self.GM actually contains some nodes with symengine.Symbol('T') instead of pycalphad.variables.T, so the result of the automatic differentiation is incorrect.

(These issues with early evaluation of properties were part of the motivation for developing the upcoming Computable Property Framework ( #432 ), which should help obviate some of these issues. This current approach will continue to work as well, though we will recommend the CPF moving forward.)

Can you try modifying your custom class like this, to see if it works?

from symengine import Symbol

class MyModel(Model):
    def build_phase(self, dbe):
        super(MyModel, self).build_phase(dbe)
        symbols = {Symbol(s): val for s, val in dbf.symbols.items()}
        for name, value in self.models.items():
            # XXX: xreplace hack because SymEngine seems to let Symbols slip in somehow
            self.models[name] = self.symbol_replace(value, symbols).xreplace(v.supported_variables_in_databases)
        self.entropyval = self.build_entropy(dbe)

    def build_entropy(self, _):
        Gibbs_energy = self.GM
        entropy = -Gibbs_energy.diff(v.T)
        return entropy

The fix for this specific issue should be to move that symbol replacement code into Model.build_phase, so that when you call the build_phase parent class, symbol replacement will occur before control is returned to your custom method. Interested in contributing a PR?

@ianhbell
Copy link
Author

ianhbell commented Jan 5, 2023

Thanks for the super fast response. Alas I must not have a new enough version since this import fails: from pycalphad.io.tdb import get_supported_variables. Do I need to build/install from source?

@richardotis
Copy link
Collaborator

Sorry about that, I was looking at a development branch. Try the edited version.

@ianhbell
Copy link
Author

ianhbell commented Jan 5, 2023

Perfecto. That hack works.

@ianhbell
Copy link
Author

ianhbell commented Jan 5, 2023

For future readers, the fixed code:

from tinydb import where
import pycalphad
from pycalphad import Model, variables as v
from pycalphad import Database
from pycalphad import calculate
from symengine import Symbol

import pandas
import numpy as np
print('Version:', pycalphad.__version__)

with open('unary50edt.tdb') as fp:
    dbf = Database(fp.read().upper())

class MyModel(Model):
    def build_phase(self, dbe):
        super(MyModel, self).build_phase(dbe)
        symbols = {Symbol(s): val for s, val in dbf.symbols.items()}
        for name, value in self.models.items():
            # XXX: xreplace hack because SymEngine seems to let Symbols slip in somehow
            self.models[name] = self.symbol_replace(value, symbols).xreplace(v.supported_variables_in_databases)
        
        self.entropyval = self.build_entropy(dbe)

    def build_entropy(self, _):
        Gibbs_energy = self.GM
        entropy = -Gibbs_energy.diff(v.T)
        return entropy

for atom in ['Ni']:
    mod = MyModel(dbf, [atom], 'LIQUID')

    # NOTICE: we need to tell pycalphad about our model for this phase
    models = {'LIQUID': mod}

    o = []
    for temp in np.arange(950, 2000, 200):
        res = calculate(dbf, [atom], 'LIQUID', P=101325, T=temp, model=models, output='SM')
        sA = res.SM.values.squeeze()
        res = calculate(dbf, [atom], 'LIQUID', P=101325, T=temp, model=models, output='entropyval')
        sB = res.entropyval.values.squeeze()
        print(temp, sA, sB)
        
    pandas.DataFrame(o).to_csv(f'{atom}_CALPHAD_test.csv', index=False)

yields

Version: 0.10.2.dev14+g8d16d69
950 74.3563867319028 74.3563867319028
1150 80.55645325671628 80.55645325671628
1350 86.13576559192114 86.13576559192114
1550 91.33373041483597 91.33373041483597
1750 96.34571598414678 96.34571598414678
1950 101.00972148214086 101.00972148214086

@ianhbell
Copy link
Author

ianhbell commented Jan 5, 2023

Is the PR simply a copy-paste into the base class in build_phase? I could try a PR, but I don't know anything about the guts.

@richardotis
Copy link
Collaborator

Yes, it's close to a copy-paste. You'd want to change the build_phase method in Model to something like:

    def build_phase(self, dbe):
        """
        Generate the symbolic form of all the contributions to this phase.
        Parameters
        ----------
        dbe : Database
        """
        contrib_vals = list(OrderedDict(self.__class__.contributions).values())
        if 'atomic_ordering_energy' in contrib_vals:
            if contrib_vals.index('atomic_ordering_energy') != (len(contrib_vals) - 1):
                # Check for a common mistake in custom models
                # Users that need to override this behavior should override build_phase
                raise ValueError('\'atomic_ordering_energy\' must be the final contribution')
        self.models.clear()
        for key, value in self.__class__.contributions:
            self.models[key] = S(getattr(self, value)(dbe))

        symbols = {Symbol(s): val for s, val in dbe.symbols.items()}

        for name, value in self.models.items():
            # XXX: xreplace hack because SymEngine seems to let Symbols slip in somehow
            self.models[name] = self.symbol_replace(value, symbols).xreplace(v.supported_variables_in_databases)

And then removing that same symbol replacement loop in the __init__ method.

ianhbell added a commit to ianhbell/pycalphad that referenced this issue Jan 5, 2023
@ianhbell ianhbell linked a pull request Jan 5, 2023 that will close this issue
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants