pycalphad
write_tdb function outputs incorrect tdb due to symengine
The write_tdb function produces an unreadable tdb due to how symengine evaluates EXP and LN. symengine.sympify is used here and is where the evaluation takes places: https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L62
I have attached a two line script which opens a tdb then writes it out again, as you can see the input (which pycalphad reads and plot correctly ) is very different to the output tdb
According to Brandon the issue may be related to this: https://github.com/symengine/symengine/pull/1828
Does pycalphad produce the same phase diagram for both TDBs?
For example this input:
FUNCTION CAO0 1 (1 - EXP(-314.8272896/T)); 10000 N !
FUNCTION CAO1 1 (1 - EXP(-125.12812657/T)); 10000 N !
FUNCTION CAO2 1 (1 - EXP(-598.37301449/T)); 10000 N !
FUNCTION CAOL0 1 (1 - EXP(-629.27558/T)); 10000 N !
FUNCTION CAOL1 1 (1 - EXP(-234.0578/T)); 10000 N !
and this output:
FUNCTION CAO0 1.0 1.0-1.0 * 2.71828182845905**(-314.8272896 * T**(-1.0));
10000.0 N !
FUNCTION CAO1 1.0 1.0-1.0 * 2.71828182845905**(-125.12812657 * T**(-1.0));
10000.0 N !
FUNCTION CAO2 1.0 1.0-1.0 * 2.71828182845905**(-598.37301449 * T**(-1.0));
10000.0 N !
FUNCTION CAOL0 1.0 1.0-1.0 * 2.71828182845905**(-629.27558 * T**(-1.0));
10000.0 N !
FUNCTION CAOL1 1.0 1.0-1.0 * 2.71828182845905**(-234.0578 * T**(-1.0));
10000.0 N !
These two should be identical for many significant figures, and I would expect pycalphad to produce identical phase diagrams for both. This is still a bug because most other Calphad implementations will not accept this output, but I want to understand if this is only a compatibility bug or if it's also a correctness bug in pycalphad.
Related: https://github.com/symengine/symengine/pull/1828#issuecomment-907627485, in particular:
We evaluate things like sin(2.0). Reason is that 2.0 is only accurate to 53 bits, so it doesn't make sense to evaluate sin(2.0) to more precision than the precision of 2.0. That's one reason to do it. The other is that it prevents explosion of the symbolic tree. So, the policy that symengine uses is that if you send in floating point number to a function like (sin, exp, etc) you get back a floating point number.
The binplot in Pycalphad doesn’t produce a plot for the output tdb but does produce one for the input.
It’s specifically the EXP causing the issue, removing those functions will produce a TDB that will plot in pycalphad.
We have logic in the TDB writer that is supposed to detect this case: https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L469-L478
We just need to work out why it isn't being triggered here.
I believe that part is working as intended, but from my small investigation its the reading of the TDB that might be causing the issue specifically the tdb_grammar() function: https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L207
which then points to: https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L107-L116
and then points to the sympify evaluation: https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L62
The logic you referenced doesn't work in this case because the expr has already been evaluated upon creation of a Database object, so the 2.71** is written out according to the else:
https://github.com/pycalphad/pycalphad/blob/ebcfbdb4dadfcce98a40db41c444a19842d8849e/pycalphad/io/tdb.py#L473-L478
It doesn't solve the general problem, but can we tweak the writing logic for Pow to treat a base number of ~2.71 as the exponential function, for writing purposes?