Problems with the smart equilibrium and kinetic solver

[taj004]

Description

I am experiencing problems with the SmartEquilibriumSolver and the SmartKineticsSolver, given transported elements, in that the pressure, temperature and chemical species blow up after a few steps. For example, in the MWE below, in the initial step the SmartEquilibriumSolver converge, but on the next step, all the species, temperature and pressure become nan.

Example

import reaktoro as rt
import numpy as np
import matplotlib.pyplot as plt

# For use random numbers later on
np.random.seed(1)

# Database
db = rt.SupcrtDatabase("supcrtbl")

# System
system = rt.ChemicalSystem(
    # The standard business
    db,
    rt.AqueousPhase().set( rt.ActivityModelHKF() ),     
    rt.MineralPhases("Calcite Magnesite Dolomite"),
    )

#%% Prepare to solve

day = 1 * 60 # one day, in seconds
duration = 5 * day # final time
steps = 30 # number of time steps
dt = duration / steps

table = rt.Table() # table to store simulation information

state = rt.ChemicalState(system)
state.set("H2O(aq)",   1.0, "kg")
state.set("Calcite",   1.0, "mol")
state.set("Magnesite", 1.0, "mol")
state.set("O2(aq)",    1.0, "umol") 

state.scalePhaseVolume("Calcite", 0.1, "cm3")

new_state = state.clone()

# Solver
solver = rt.SmartEquilibriumSolver(system)
cond = rt.EquilibriumConditions(system)

options = rt.SmartEquilibriumOptions()
#options.learning.epsilon = 1e-7
#options.learning.optima.output.active=True
options.learning.optima.convergence.tolerance = 1e-6

solver.setOptions(options)

# Initial elements
b = new_state.componentAmounts().asarray()
components = state.componentAmounts().asarray().size 
#print(f"Initial elements {b}")

#%% Solve

# Time loop
for i in range(steps + 1):
    print(i)
    # Let this represent `transported' elements
    b +=  5e-7 * np.random.uniform(components)
   
    cond.setInitialComponentAmounts(b)
    
    # Solve the current system
    result = solver.solve(new_state, cond) 
    
    assert result.succeeded(), f"calculation faild at time step {i}" 
    
    # Get current information about the system
    props = new_state.props() 
    aprops = rt.AqueousProps.compute(props)
    
    # Store in table
    table.column("Time") << i*dt / day # days
    
    table.column("Calcite") << props.speciesAmount("Calcite").val()
    table.column("Magnesite") << props.speciesAmount("Magnesite").val()
    table.column("Dolomite") << props.speciesAmount("Dolomite").val()
    table.column("Ca+2") << props.speciesAmount("Ca+2")
    table.column("Mg+2") << props.speciesAmount("Mg+2").val()
       
    table.column("ph") << aprops.pH() 
    
    table.column("H") << props.elementAmount("H").val()
    table.column("C") << props.elementAmount("C").val()
    
# end time-loop     

#print(f"Elements at last step {b}")

# Compare states
# print("Initial state")
# print(state)
# print("Final state")
# print(new_state)

#%% Plot

plt.figure()
plt.plot(table["Time"], table["Ca+2"], color="blue", label="$Ca^{2+}$")
plt.plot(table["Time"], table["Mg+2"], color="orange", label="$Mg^{2+}$")
plt.ylabel("Amount [mol]")
plt.xlabel("Time [day]")
plt.legend(loc="upper right")
plt.show()   
 
plt.Figure()
plt.plot(table["Time"], table["H"], label="H")
plt.ylabel("Amount [mol]")
plt.xlabel("Time [day]")
plt.legend(loc="upper right")
plt.legend(loc="upper right")

What is your operating system?

Linux

Additional Information

No response

[allanleal]

Hi @taj004 , thanks for reporting this issue and providing an example. I'll check/fix this issue next week.

[allanleal]

Hi @taj004 , could you please show me the versions of the packages installed in your conda environment? In the terminal, with the conda environment containing reaktoro activated, please do:

conda list

On my end, I had to set temperature and pressure in the EquilibriumConditions object (otherwise runtime errors were raised telling me that temperature and pressure were not set). Whenever EquilibriumConditions are used, and temperature and pressure are know in the problem, they must be specified.

cond.temperature(25.0, "°C")
cond.pressure(1.0, "bar")
cond.setInitialComponentAmounts(b)

The figures produced were:

SmartEquilibriumSolver and SmartKineticsSolver are still experimental. Feel free to report any issues!

[taj004]

Hello @allanleal. Thanks for the answer. Since my reaktoro conda environment has numerous packages, I only show those related to reaktoro (autodiff, optima etc; If have not included an important one, please let me know.). I use a bit older version since I experienced the same issues as discussed in #344 with the newer ones. Anyways, your solution of setting temperature and pressure in cond solved the problem

[allanleal]

Hi @taj004 , thanks for providing the information above. Newer versions of Reaktoro raise an error if temperature and pressure are not specified in EquilibriumConditions when they are imposed constraints. Older versions cause this obscure error you were facing (and reported here).

I understand that in the current version of Reaktoro (v2.9.3) for some calculations the configuration:

options EquilibriumOptions()
options.optima.backtracksearch.apply_min_max_fix_and_accept = True

solver = EquilibriumSolver(specs)
solver.setOptions(options)

is necessary to recover the previous algorithm behavior (otherwise the calculation fails), as reported in #344 .