pH Calc slow inside PDEs solver

[pycckuu]

Hi,

I would like to use your “pH Calc” algorithm in my Reactive-Transport-Toolbox.

This reaction toolbox solves advection-diffusion-reaction PDE for saturated (at the moment) and not saturated (will be added soon) flow. The problem that I encounter at the moment that I solve the pH for vectors and I am iterating across indexes of concentrations. In this implementation of pHcalc it takes very long time to solve for each time step in PDE. Thus, I was wondering if it is possible to accelerate the estimation of pH in your toolbox? One of the solution would be to provide a jacobian. Would it be possible to make?

Thanks.

Looking forward to your reply.

My toolbox is here: https://github.com/biogeochemistry/PorousMediaLab

[rnelsonchem]

@pycckuu Thanks for the interest in the project. Please see this Gist for some suggestions for speeding up the calculations. I hope that one of these helps you out.

Unfortunately, I won't have time myself to modify the code to allow the use of a jacobian for the solver, if that is really what you need. But I think if you take a look at the code, it probably shouldn't be too difficult to implement. The scipy.optimize.minimize function does accept a jacobian already, so you could just modify System._diff_pos_neg as a start.

All the best

[pycckuu]

Thanks for the gist, @rnelsonchem

Yes, I tried to input the jacobian and I also used the numerical approximation of one scipy.optimize.approx_fprime. Both of them did not really help to improve the computation time. Moreover, they were not quite stable and were very dependent on the initial guess. Thus, the simplex method of Nelder-Mead still was the most robust one.

The suggestion to use a narrower pH window is really great given that for the pH values we are not require the precision greater that 0.01. Moreover, I think we can use numpy.vectorize here to speed up the calculations on vectors a bit. I will try it soon and let you know how it goes.

By the way, I was wondering how difficult would it be to include the buffering/neutralizing systems in this approach? For example, CaCO3(s) base with corresponding pKs which can precipitate or dissolve based on pH (see attached picture as an example). I was thinking about this and could not find an elegant way (in a manner how you did with the acids) to incorporate it in the system. Do you have any ideas? I can implement it if you have.

Thanks

[rnelsonchem]

In principle, this methodology would be appropriate for the system you describe. However, I coded this module specifically for simple aqueous acid/base systems, so you wouldn't be able to directly use this code for that problem. The problem is that the Ca2+ concentration (a Neutral in pHcalc) is dependent on the concentration of one component from a polyprotic acid (an Acid in pHcalc). I've thought about this, and I think you would need to essentially rewrite the entire module from scratch. It would also be nice to have an estimate of ionic strength as well. One of these days, perhaps.