Latest reaktoro version 2.9.3 shows warnings and is slower compared to previous calculations

[bgvinayk]

I have been running equilibrium calculations in latest version as shown in image below. It runs slower compared to previous version and gives warnings.
I tried to install older version but getting vector too long error.

Please help regarding this issue. Thank you.

1. Is there a way to run in latest version faster
2. Is there a way to handle vector too long error
3. Is there a proper way to downgrade the version without any dependency issues

from reaktoro import *
import warnings
warnings.filterwarnings('ignore')

db = PhreeqcDatabase("minteq.v4.dat")

solution = AqueousPhase(speciate("H O Ca Mg Na K Zn Fe Mn C Si P S Cl F N Ba Sr Al Cu Mo"))
solution.set(ActivityModelDebyeHuckelPHREEQC())

# gases = GaseousPhase("CO2(g)")
# gases.setActivityModel(ActivityModelPengRobinson())

system = ChemicalSystem(db, solution)

specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()

solver = EquilibriumSolver(specs)
state = ChemicalState(system)

Edit on 30.10.23 by @allanleal: Using syntax highlighting for given Python example.

[allanleal]

Hi @bgvinayk, thanks for reporting this issue.

Chemical equilibrium calculations can fail for a variety of reasons (for example, input conditions may be outside the range of validity of some internal thermodynamic models, or the algorithm has difficulty converging with certain initial conditions).

You provide an example above, but it does not contain the details necessary to reproduce the failure. Could you please review this and provide a complete example of a failed chemical equilibrium calculation?

1. Is there a way to run the latest version faster

What may be happening is that some changes in the algorithm is causing errors in your calculations that were not happening before. The failed calculation continues up to the maximum number of iterations (currently 200 by default). You then get the impression that the latest version is not as fast as before. The example requested above would be very helpful to find out.

2. Is there a way to handle very long vector errors

I suspect that you installed an older version of Reaktoro, but still rely on a newer version of Optima in the same environment. This may be the cause of the vector too long problem. You would need to downgrade Optima too, but that should have happened automatically. What is the older version of Reaktoro you have installed?

3. Is there a proper way to downgrade the version without dependency issues

See above.

I look forward to receiving a complete example of your failing equilibrium calculation. If you want to reverse the algorithmic behavior of Reaktoro v2.9 to that of v2.8 and before, do the following:

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

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

[bgvinayk]

Thanks Allan, your solution worked. For your reference, these are my inputs. These worked in my previous 2.7.3 version. Not sure why there are outside thermodynamic conditions.

from reaktoro import *
import warnings
warnings.filterwarnings('ignore')

db = PhreeqcDatabase("minteq.v4.dat")

solution = AqueousPhase(speciate("H O Ca Mg Na K Zn Fe Mn C Si P S Cl F N Ba Sr Al Cu Mo"))
solution.set(ActivityModelDebyeHuckelPHREEQC())

# gases = GaseousPhase("CO2(g)")
# gases.setActivityModel(ActivityModelPengRobinson())

system = ChemicalSystem(db, solution)

specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()

solver = EquilibriumSolver(specs)
state = ChemicalState(system)

state.set("H2O", 1.0, "kg")
state.set("Na+", 1.0, "mol")
state.set("Cl-", 1.0, "mol")
# state.set("CO2", 0.4, "mol")
state.set("Cu+2", 0.01, "mol")
state.set("Ba+2", 0.000001, "mol")
state.set("Sr+2", 0.000001, "mol")
state.set("SO4-2", 0.001, "mol")

conditions = EquilibriumConditions(specs)
conditions.temperature(40.0, "celsius")
conditions.pressure(1, "atm")
conditions.pH(8.5)
solver.solve(state, conditions)

aprops = AqueousProps(state)

[allanleal]

Hi @bgvinayk , good to know that changing that algorithm option resolved your issue. I'll be able to investigate the numerical issue on the calculation above (using Reaktoro v2.9+) in the next days.

"Not sure why there are outside thermodynamic conditions."

Most likely the issue is numerical. However, I've noticed in the example above that charge is far from balanced. Although Reaktoro can handle non-zero charge solutions, in some cases significant charge imbalance can result in numerical convergence difficulties. Setting both SO4-2 and Cu+2 to 0.001 mol would significantly improve charge balance. Not sure if your original compositional values are intentional.

[allanleal]

On a separate issue, instead of ActivityModelDebyeHuckelPHREEQC, you may want to use ActivityModelPhreeqc, which is consistent with how PHREEQC evaluates activity coefficients of water and solutes when using databases such as phreeqc.dat, llnl.dat, minteq.v4.dat. For pitzer.dat, use ActivityModelPitzer.

[bgvinayk]

I tried to use ActivityModelPhreeqc instead of ActivityModelDebyeHuckelPHREEQC, but was facing an error

solution.setActivityModel(ActivityModelPhreeqc())
TypeError: ActivityModelPhreeqc(): incompatible function arguments. The following argument types are supported:
1. (arg0: reaktoro.reaktoro4py.PhreeqcDatabase) -> Callable[[reaktoro.reaktoro4py.SpeciesList], reaktoro.reaktoro4py.ActivityModel]

Invoked with:

[allanleal]

You need to pass the Database object to the function: ActivityModelPhreeqc(db)

[bgvinayk]

Hi @allanleal,
I am attaching a sample code below to indicate the issue. It takes more time to run pH=8.8 scenario with warnings. But scenarios of pH = 8.7 and pH=8.9 run fast and without warnings. Please suggest how I can improve performance for a generic scenario. Appreciate your prompt response. Thanks!

from reaktoro import *
import time

startTime = time.time()

db = PhreeqcDatabase("minteq.v4.dat")

solution = AqueousPhase(speciate("H O Ca Mg Na K Zn Fe Mn C Si P S Cl F N Ba Sr Al Cu"))
solution.setActivityModel(ActivityModelDebyeHuckelPHREEQC())

gases = GaseousPhase("CO2(g)")
gases.setActivityModel(ActivityModelPengRobinson())

system = ChemicalSystem(db, solution, gases)

specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()

options = EquilibriumOptions()
options.optima.backtracksearch.apply_min_max_fix_and_accept = True
solver = EquilibriumSolver(specs)
solver.setOptions(options)

water_mass = 1000
P = 1
Tc = 40
'pH = 8.7'
pH = 8.8
'pH = 8.9'

state = ChemicalState(system)
state.pressure(P, "atm")
state.temperature(Tc, "celsius")
state.set("H2O"    , water_mass, "kg")
state.set("Na+"    , 2, "mol")
state.set("K+"    , 0, "mol")
state.set("Ca+2"    , 54, "g")
state.set("Mg+2"    , 9, "g")
state.set("Ba+2"    , 0, "g")
state.set("Sr+2"    , 0.2, "g")
state.set("Fe+2"    , 0, "g")
state.set("Zn+2"    , 0, "g")
state.set("Al+3"    , 0, "g")
state.set("Mn+2"    , 0, "g")
state.set("Cu+2"    , 0, "g")
state.set("NH3"    , 0, "g")
state.set("CO3-2"    ,  2.5, "mol")
state.set("H4SiO4"    ,  1.3, "mol")
state.set("SO4-2"   , 1, "mol")
state.set("PO4-3"   , 0, "g")
state.set("Cl-"    , 1.34, "mol")
state.set("NO3-"    , 0, "g")
state.set("F-"    , 0, "g")

conditions = EquilibriumConditions(specs)
conditions.temperature(Tc, "celsius")
conditions.pressure(P, "atm")
conditions.pH(pH)

result = solver.solve(state, conditions)

print("Duration: ", time.time()-startTime, "seconds")

[bgvinayk]

@allanleal I ran the same code and got following in the output. Not sure why am I getting the warning and higher time. Can you suggest and what version to use.

[allanleal]

Hi @bgvinayk, I used Reaktoro v2.9.2, Linux. The longer time is not due to the warning message, but because the calculation took 200 iterations before stopping (as unsuccessful). Let me know the version and OS you used please.

[bgvinayk]

@allanleal Version 2.9.3 and windows 10. Case pH 8.7 and pH 8.9 has no warning and 5 times less time compared to pH 8.8 case.

[allanleal]

Hi @bgvinayk, I ran this on Windows, v2.9.4. The calculation fails for pH = 8.8 and no warning is given (which should). I'll need to investigate further. There are some planned improvements to the algorithm that would solve this issue. I'll be happy to explain more in a short call.

[bgvinayk]

Thanks @allanleal for your assistance on this. I understood your explanation. Will look forward to improvements in coming releases.