Example of a problem in which the scale of amounts cause the equilibrium algorithm to fail

[allanleal]

The example below works when rate = 1 but fails when rate = 1000:

from reaktoro import *

db = NasaDatabase("nasa-cea")

solid = SolidPhase("C(gr)")
condensed = CondensedPhase("Fe.947O(cd) Fe(cd)")
gas = GaseousPhase("O2 N2 H2O CO CO2 H2")

system = ChemicalSystem(db, gas, solid, condensed)

rate = 1      # THIS WORKS
# rate = 1000   # DOES NOT WORK
T = 300

state = ChemicalState(system)
state.set('Fe.947O(cd)', rate, 'mol')
state.set('C(gr)', rate, 'mol')
state.temperature(T, 'celsius')
state.pressure(5, 'atm')
solver = EquilibriumSolver(system)

options = EquilibriumOptions()
options.optima.output.active = True
solver.setOptions(options)

result = solver.solve(state)

assert result.succeeded()

The attached files below show the progress of the equilibrium calculation. The case rate = 1000 should (ideally) behave the same as the case rate = 1.

optima.log-rate-1.txt optima.log-rate-1000.txt

TIP: In the meantime, while a fix for this issue is produced, prefer working with normalized input conditions (e.g., normalize initial amounts of species by a total amount) and then scale back the calculated species amounts by the normalizing factor.

Thanks, @Noxical, for raising this issue via Gitter.

[allanleal]

This issue still exist, despite some recent improvements in the algorithm due to scale of the variables. I've noticed that changing the lower bounds from 1e-16 to 1e-10 for this problem, with high species amounts, caused the calculation to converge without issues. Tested with Reaktoro v2.8.2.

from reaktoro import *

db = NasaDatabase("nasa-cea")

solid = SolidPhase("C(gr)")
condensed = CondensedPhase("Fe.947O(cd) Fe(cd)")
gas = GaseousPhase("O2 N2 H2O CO CO2 H2")

system = ChemicalSystem(db, gas, solid, condensed)

# rate = 1      # THIS WORKS
rate = 1000   # THIS WORKS WITH LOWER BOUNDS FOR SPECIES AMOUNTS AS 1E-10
T = 300

state = ChemicalState(system)
state.set('Fe.947O(cd)', rate, 'mol')
state.set('C(gr)', rate, 'mol')
state.temperature(T, 'celsius')
state.pressure(5, 'atm')
solver = EquilibriumSolver(system)

options = EquilibriumOptions()
options.optima.output.active = True
options.epsilon = 1e-10  # set lower bounds to 1e-10 moles instead of 1e-16
solver.setOptions(options)

result = solver.solve(state)

assert result.succeeded()