Terminate IPhreeqc.dll from matlab during 1D transport

[Mouadh]

Dear Phreeqc users,


I'm a postdoc working on reactive transport modelling of CO2 interaction with basalt.


I'm running 1D transport in phreeqc using the IPhreeqc.dll library from matlab.
 My model works fine, but I 'm trying to optimize some parameters in which some time phreeqc fails to converge to a solution.
 The issue is that I have more than 1000 shifts in my transport process and the model takes forever if the tested parameters gives convergence error.
Is it possible to add a time limit for the transport runtime?
 
I tried to find a way to do this from matlab but according to this post it is not really possible (https://se.mathworks.com/matlabcentral/answers/454098-how-to-do-wait-calllib-to-give-answer-or-catch-exception)
Thank you!
Mouadh Addassi

[dlparkhurst]

There is no feature in PHREEQC to end a calculation based on clock time.

I don't completely understand what is causing your problem. On the one hand, you have 1000 shifts, which could take a while, but it would take a while whether the calculation succeeds or fails. It only takes failure in one cell for the entire calculation to fail, and relatively, it would not take too much longer to have one calculation fail than to do the entire calculation.

Options in the KNOBS data block may or may not help with your problem. I generally use the default parameters. The only options I tend to change are -step (say 10) and -pe_step (say 5), which may be more likely to cause convergence. If you make the values too small, you need more -iterations, and then it presumably will take even longer if it continues to fail. I don't know your system, but very small concentrations can be a problem, like iron in equilibrium with a ferric oxyhydroxide mineral. In this case, -diagonal_scale sometimes helps. One other common problem occurs in solutions without redox buffering, in this case, adding a redox species can help: SOLUTION_SPECIES; H2O + 0.01e- = H2O-0.01; log_k -8

An undocumented option is -try n, where n is between 1 and about 15. PHREEQC has a set of convergence parameters that are varied if the initial set fails to produce convergence. In all there are about 15 sets of parameters that are tried. So when convergence fails completely, the calculation has been attempted with all sets of parameters; further, the number of iterations is doubled after the first set, which also means that the calculation takes longer. So, if you use -try 4, only 4 sets of convergence parameters are tried, and the calculation should fail more quickly. But again, I don't see how a failed calculation in one cell would cause the entire calculation to be slow. It should be no worse than calculating about 50 cells, which would be small relative to your entire calculation.

If you can attach a simple input file (and database if nonstandard) that demonstrates the problem, I may be able to give a better evaluation of your problem.

[Mouadh]

Thank you for a quick response. I will try adding redox buffering. I 'm already trying different sets of  convergence parameters (matlab file attached)

My problem is that the calculation does not fail completely but it becomes much much slower. Normally one iteration takes less than 30 seconds (for example if log_k basalt_gl = 1), but in cases like the attached one (log_k basalt_gl = -1), it takes more than 20 minutes.

Al  1e-10 in solution 1-7 was added as an attempt to help the problem but it did not change much.

[dlparkhurst]

I don't think it is really a convergence problem. Try adding some Al and Si (say 1e-10) to SOLUTION 0 and all of the initial cell solutions. Without Al and Si, it is not possible to calculate SR for the basalt, and I think it creates a discontinuity that is hard to integrate in KINETICS. The -cvode option fails entirely, but runs successfully after some Al and Si is added. Hopefully, you can simply use the default convergence parameters.

You can test whether you want to use the default Runge Kutta method or CVODE for KINETICS. For the file that you sent, cvode ran slightly faster, but usually Runge Kutta is faster,. If you start to add more kinetic reactions (creating a "stiff" system of rates), cvode may be the better option.

[Mouadh]

 Thank you for taking time out of your personal time to help me.  I am truly grateful.

Adding Al and Si to all Solutions helps a lot.

Best,
Mouadh