Root finding failed
zplizzi opened this issue · 4 comments
When trying to find the equilibria of the following system, I get "root finding failed" and "too much of at least one component" warnings.
eqsys = EqSystem.from_string("""
CO2 + H2O = H2CO3; 1.7e-3
H2CO3 = H+ + HCO3-; 10**-6.3
HCO3- = H+ + CO3-2; 10**-10.3
HPO4-2 + H2O = H2PO4- + OH-; 10**-6.79
Na2HPO4 = 2 Na+ + HPO4-2; 10**5
KH2PO4 = K+ + H2PO4-; 10**5
HCl = H+ + Cl-; 10**5.9
CH3CO2H = CH3CO2- + H+; 10**-4.756
H2O = H+ + OH-; 10**-14
""")
arr, info, sane = eqsys.root(defaultdict(float, \
{'H2O': 55.4,
'CO2': 3,
"KH2PO4": 0.0422 * 1,
"Na2HPO4": 0.0578 * 1,
"HCl": 0,
"CH3CO2H": 0.0,
"CH3CO2-": 0
}))
conc = dict(zip(eqsys.substances, arr))
print(conc)Even with somewhat simplified versions of this system, I still get the errors at times.
Hi, it was a while since I was solving equilibria, but I remember that the solvers provided by scipy had trouble with convergence quite often. I tried applying various transformations with some success, but I think using another solver is probably a better solution long term, I started looking into using IPOPT (cyipopt provides Python bindings), but I got busy with other work.
For trying variable transformations, I have an example notebook (ammonical cupric solution), if you want to try IPOPT there will be some development that needs to be done, I won't have time to do this myself in the foreseeable future, but if you're interested in trying that avenue, let me know and I'll try to help out.
Either way, you'll probably want to indicate that some of these species belong to another phase (it's solubility products for the potassium and sodium salts I presume?). I seem to remember that the default parser respected the (s) suffix?
No worries, I probably don't have the time/skill to dig into improving the solver. I'll just work with what you have here - this library is awesome overall!
I'm happy to hear that you are finding the library useful! And I appreciate that you're opening issues as well, it will make it easier to prioritize what to improve when I find some spare time.
Just saw that openturns wraps Ipopt in their python wrapper, leaving a link here for future reference:
http://openturns.github.io/openturns/latest/user_manual/_generated/openturns.Ipopt.html#ipopt