VCS equilibrium solver jumps between different solutions for mixtures with binary-solution-tabulated phase

[spinnau]

Problem description

Calling the equilibrate method on mixtures containing a `binary-solution-tabulated phase using the VCS solver can give strange results. I have noticed that there are jumps between different solutions when the solver is called several times.

This is not a general problem of the VCS solver. It works fine even with other more complex models e.g. containing HMW-electrolyte phases. When changing the binary-solution-tabulated phase to ideal-condensed, this problem doesn't occur.

Steps to reproduce

yaml = """
phases:
- name: Li(s)
  thermo: fixed-stoichiometry
  elements: [Li]
  species: [Li(cr)]
  state: {T: 298.15, P: 1.0e+5 Pa}

species:
- name: Li(cr)
  composition: {Li: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 453.69]
    data:
    - [0.610909942, 0.0141041217, -1.7495817e-05, -3.33741023e-08, 7.76629665e-11,
      -625.121208, -3.26449947]
    note: TPIS82
"""

p1 = ct.Solution('lithium_ion_battery.yaml', 'anode')
p2 = ct.Solution(yaml=yaml, name='Li(s)')
mix = ct.Mixture([p1, p2])

mix.T = 293
mix.P = 1e5
mix.species_moles = [1.0, 1.0, 0.0]

mix.equilibrate('TP', solver='vcs', log_level=1)
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

Behavior

The successive calls to the equilibrate method produce different results. Without changing Temperature, pressure or composition the results should be the same:

Trying VCS equilibrium solver
VCS solver succeeded
[0.04620062 1.95379938 0.95379938]
[1. 1. 0.]
[0.04620062 1.95379938 0.95379938]
[1. 1. 0.]
[0.04620062 1.95379938 0.95379938]

I've also changed the rtol and estimate_equil parameters, but the problem remains the same.

System information

• Cantera version: 3.1.0
• OS: Arch Linux
• Python 3.13.5

Additional context

The solution of the VCS solver is also different from the Gibbs solver. And for this example, the VCS solver gives the same result if the equilibrate method is called again after using the Gibbs solver:

mix.equilibrate('TP', solver='gibbs')
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

mix.equilibrate('TP', solver='vcs')
print(mix.species_moles)

[0.2347416 1.7652584 0.7652584]
[0.2347416 1.7652584 0.7652584]
[0.2347416 1.7652584 0.7652584]

In principle this looks like I could just use the Gibbs solver, but that very often gives errors like in #1023.

[speth]

Thanks for reporting this. I'd add that it's not clear which if any of these "solutions" is actually correct. If you print out the total chemical potential for the mixture:

for _ in range(5):
  mix.equilibrate('TP', solver='vcs')
  print(f'VCS: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

mix.equilibrate('TP', solver='gibbs')
print(f'Gibbs: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

You will get the following output:

VCS: n = [0.04620062 1.95379938 0.95379938], G = -1.8996e+07
VCS: n = [1. 1. 0.], G = -1.1984e+07
VCS: n = [0.04620062 1.95379938 0.95379938], G = -1.8996e+07
VCS: n = [1. 1. 0.], G = -1.1984e+07
VCS: n = [0.04620062 1.95379938 0.95379938], G = -1.8996e+07
Gibbs: n = [0.2347416 1.7652584 0.7652584], G = -1.8337e+07

Which shows that the lowest value is the first one returned by the VCS solver, while the solution returned by the Gibbs solver has a higher chemical potential, so is clearly not correct.

To investigate this a bit further, we can use the fact that this example conveniently has a single degree of freedom, which can be characterized as the number of moles of Li[Anode], with the other two species moles being fully constrained by the element conservation constraints. If we plot the chemical potential as a function of this one free variable, we can see what the space looks like:

import cantera as ct
import numpy as np
import matplotlib.pyplot as plt

yaml = """
phases:
- name: Li(s)
  thermo: fixed-stoichiometry
  elements: [Li]
  species: [Li(cr)]
  state: {T: 298.15, P: 1.0e+5 Pa}

species:
- name: Li(cr)
  composition: {Li: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 453.69]
    data:
    - [0.610909942, 0.0141041217, -1.7495817e-05, -3.33741023e-08, 7.76629665e-11,
      -625.121208, -3.26449947]
    note: TPIS82
"""

p1 = ct.Solution('lithium_ion_battery.yaml', 'anode')
p2 = ct.Solution(yaml=yaml, name='Li(s)')
mix = ct.Mixture([p1, p2])

mix.T = 293
mix.P = 1e5
mix.species_moles = [1.0, 1.0, 0.0]

for _ in range(5):
  mix.equilibrate('TP', solver='vcs')
  print(f'VCS: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

mix.equilibrate('TP', solver='gibbs')
print(f'Gibbs: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

x = np.linspace(0,1,2000)
gg = np.zeros(x.shape)

for i in range(len(x)):
    mix.species_moles = [x[i], 2-x[i], 1-x[i]]
    gg[i] = mix.species_moles @ mix.chemical_potentials

fig, ax = plt.subplots()

ax.plot(x, gg)

mix.species_moles = [1.0, 1.0, 0.0]
for _ in range(2):
    mix.equilibrate('TP', solver='vcs')
    ax.plot(mix.species_moles[0], mix.species_moles @ mix.chemical_potentials, 'o', label='VCS')

mix.equilibrate('TP', solver='gibbs')
ax.plot(mix.species_moles[0], mix.species_moles @ mix.chemical_potentials, '*', label='gibbs')
ax.set(xlabel='Li[anode] moles', ylabel='Chemical potential [J/kmol]')
ax.legend()
fig.savefig('issue1928.png', dpi=200)
plt.show()

Which gives the following plot:

Where you can see that the solution from the gibbs solver is clearly wrong, and isn't even at a local minimum, and even the "good" solution from the VCS solver seems to be slightly off from the actual minimum of the chemical potential. The other thing clearly visible here is that the function is fairly lumpy, driven by the discrete nature of the "tabulated" thermo data. This lack of smoothness may be driving some of the problematic behavior for both solvers.

[spinnau]

Thanks for looking into this and your nice analysis of the problem. I have tried to smooth the tabulated data, but unfortunately that changes almost nothing.

import cantera as ct
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from ruamel.yaml import YAML


yaml = """
phases:
- name: Li(s)
  thermo: fixed-stoichiometry
  elements: [Li]
  species: [Li(cr)]
  state: {T: 298.15, P: 1.0e+5 Pa}

species:
- name: Li(cr)
  composition: {Li: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 453.69]
    data:
    - [0.610909942, 0.0141041217, -1.7495817e-05, -3.33741023e-08, 7.76629665e-11,
      -625.121208, -3.26449947]
    note: TPIS82
"""

p1 = ct.Solution('lithium_ion_battery.yaml', 'anode')
p2 = ct.Solution(yaml=yaml, name='Li(s)')

tab = p1.input_data['tabulated-thermo']
x_tab = tab['mole-fractions']
x = np.linspace(min(x_tab), max(x_tab), 4000)
#x = np.linspace(0, 1, 4000)


fig, axs = plt.subplots(1, 3, figsize=(10,4), sharex=True)

yaml = YAML().load(p1.write_yaml())
yaml['phases'][0]['tabulated-thermo']['mole-fractions'] = x.tolist()

for i, var in enumerate(['enthalpy', 'entropy', 'molar-volume']):
    y_tab = p1.input_data['tabulated-thermo'][var]
    interp = sp.interpolate.CubicSpline(x_tab, y_tab, extrapolate=True)
    y = interp(x)
    
    yaml['phases'][0]['tabulated-thermo'][var] = y.tolist()

    axs[i].plot(x_tab, y_tab, marker='o', label='tab')
    axs[i].plot(x, y, label='interpolated')
    axs[i].set_ylabel(var)

axs[0].legend()
fig.tight_layout()
plt.show()


p1_smooth = ct.Solution(yaml=str(yaml), name='anode')
mix = ct.Mixture([p1_smooth, p2])

mix.T = 293
mix.P = 1e5
mix.species_moles = [1.0, 1.0, 0.0]

for _ in range(5):
  mix.equilibrate('TP', solver='vcs')
  print(f'VCS: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

mix.equilibrate('TP', solver='gibbs')
print(f'Gibbs: n = {mix.species_moles}, G = {mix.species_moles @ mix.chemical_potentials:.4e}')

x = np.linspace(0,1,2000)
gg = np.zeros(x.shape)

for i in range(len(x)):
    mix.species_moles = [x[i], 2-x[i], 1-x[i]]
    gg[i] = mix.species_moles @ mix.chemical_potentials

fig2, ax = plt.subplots()

ax.plot(x, gg)

mix.species_moles = [1.0, 1.0, 0.0]
for _ in range(2):
    mix.equilibrate('TP', solver='vcs')
    ax.plot(mix.species_moles[0], mix.species_moles @ mix.chemical_potentials, 'o', label='VCS')
    for ax_ in axs: ax_.axvline(mix.species_moles[0], color='grey')


mix.equilibrate('TP', solver='gibbs')
ax.plot(mix.species_moles[0], mix.species_moles @ mix.chemical_potentials, '*', label='gibbs')
for ax_ in axs: ax_.axvline(mix.species_moles[0], color='red')
ax.set(xlabel='Li[anode] moles', ylabel='Chemical potential [J/kmol]')
ax.legend()
fig2.savefig('issue1928_2.png', dpi=200)
plt.show()