The Crashing in Matlab programming

[sharmilan]

Dear Phreeqc experts,

For my reactive transport model, I am using PhreeqcRM.dll library (phreeqcrm-3.6.2-15100) with call function in Matlab for chemical equilibrium.

when I am running chemical equilibrium in PhreeqcRM, if there is any failure in any cell, my Matlab program is crashed immediately. Therefore, I have to restart my Matlab program from initial.   

But, when I am using IPhreeqc, if there is any failure in any cell, the Matlab program does not crash and also SelectedOutputValue gives unreasonable value with an error message. So, I can able to continue my programme with the previous step value (without reinitializing Matlab program).

Therefore, is there any option or way in PhreeqcRM to avoid the crashing in Matlab programming?

Thank you
sharmilan.

[dlparkhurst]

Looks like there is a try/catch capability in Matlab.

[sharmilan]

Dear Sir

Thank you for your reply, my PhreeqcRM model crashed in after some step. the reason for failure does not print in the output file. do we need to use any key function to print the reason for failure?

here, I am attaching output file also.
Node: total cell nos: 31, but the first two nodes results are set for printing in the output file.   

Thank you
Sharmilan.S

[sharmilan]

Let me clarify the scenario more clearly.

My PhreeqcRM model exits with the error message "Unknown Exception" when I load the input parameters using 'RM_RunString' (input parameters are attached) at the first time step.
1) The reason for failure does not reflect in the output log file. However, I have already set the 'RM_SetFilePrefix' and 'RM_OpenFiles'.  Do we need to use any function to initiate the logging process?
2) Also when PhreeqcRM throws the error message, it stops my script entire script. I tried to bypass the exception using try/catch blocks. Unfortunately, Matlab does not treat it as an exception. Instead of jumping to catch block, it exits the code.

Thank you
Sharmilan

[dlparkhurst]

Try 3 cells and print all output with SetPrintChemistryOn.

I have begged you to simplify and add complexity sequentially to see where problems occur. You have chosen not to do this, or at least not relay these attempts, and so, I am not sympathetic to your problems.

I cannot help you debug this, unless you attach a PHREEQC input file with associated database that fails. A C++ program that fails would also work.

I am skeptical of your approach where you add large mole transfers through REACTION. If this is a flux between cells, I think you need to take smaller time steps so that the REACTION transfers are not so large. If you are somehow aggregating all the reactants into REACTION for a cell, I think that is not the right approach.

[dlparkhurst]

Reading through your input file, it seems that you start with a kilogram of water and add all the solutes and solids with REACTION. It is not clear to me that you conserve water, because you always start with a kilogram of water and do not adjust it in the reaction that I can see. But that is your problem.

PhreeqcRM was designed to keep the solids in a cell and transport the water and solutes. It may have trouble when 20 moles of Ca is thrown in and it has to find the stable phase assemblage among many excess minerals, surfaces, and, solid solutions.

At the nth step, it appears you start with 1 kg water and add a reaction, R(n), and then calculate equilibrium phases, surfaces, and solid solutions. You then transport and at step n+1, you again start with a kilogram of water and add R(n+1), setting equilibrium phases, surfaces, and solid solutions to 0. I think the calculations would be more stable, if you do not zero out the equilibrium phases, surfaces, and solid solutions, but rather leave them intact and at step n+1, add the reaction R(n+1) - R(n). That would be closer to the way PhreeqcRM was designed to be used. It would simplify the calculation of the stable phase assemblage and avoid the possibility of dealing with 20 molal Ca in the intermediate sequential iterations.

[sharmilan]

Dear Sir

Thank you for your reply, as you mentioned in the previous mail, I used smaller time steps for transfer mol in REACTION keyword and I found the problem in my solid solution phase (CSHQ).  when I am running my model without surface complexation, my model is working without any failure in chemical equilibrium. but when I include the surface complexation, my PhreeqcRM model crashed in after some step.

Due to the failure of chemical equilibrium model, PhreeqcRM throws the error message and it stops my entire Matlab script. I tried to bypass the exception using try/catch blocks. Unfortunately, Matlab does not treat it as an exception. Instead of jumping to catch block, it exits the code. Therefore, in PhreecRM, is there any option to only throws the error message without stopping my entire Matlab scrip?  (i know, it was possible when I use Phreeqc previously)

Thank you
Sharmilan.S

[dlparkhurst]

Are you still adding all solutes and reactants with REACTION? I do not think PhreeqcRM will handle that. You will need to use the R(n+1) - R(n) approach.

There are 3 options for error handling with the PhreeqcRM method SetErrorHandlerMode. I don't know if Matlab will be happy with any of them.

"Set the action to be taken when the reaction module encounters an error. Options are 0, return to calling program with an error return code (default); 1, throw an exception, in C++, the exception can be caught, for C and Fortran, the program will exit; or 2, attempt to exit gracefully."

[sharmilan]

Dear Sir

Thank you for your reply, for example, if the 10 mol of solute needs to be added in a particular time step, I use 100 iterations in each step to add that 10 mol. therefore, in each iteration, 0.1 mol (10/100) will be added through REACTION keyword with the same cell no (Note: I am not using REACTION MODIFY keyword).

Due to the carbonation of concrete analyses,
after a few time steps later, some of EQ and SS phase does not need in my chemical system, therefore, is there any method to remove or delete that particular phases from my chemical system at an intermediate timestep?  (because these excess phases are creating convergence problem in my PhreeqcRM)

Thank you
Sharmilan.S

[dlparkhurst]

Not sure how you are adding the reaction in 100 steps. If you should are using " in 100 steps", make sure INCREMENTAL_REACTIONS true; otherwise, the program will add 0.1 moles to the solution, start over and add 0.2 moles of reaction to the solution,..., and finally use the solution and add 10 moles.

I reiterate, that I think it would be best to keep the state from the last time step. The moles of solids and other reactants will be nearly correct (and changes would diminish with smaller time steps); it will only be necessary for PHREEQC to calculate adjustments, rather than starting from pure water and having to calculate the stable phase assemblage from scratch. Calculating the difference between two REACTION definitions is relatively easy. You can append the previous REACTION definition (without "REACTION") to the current REACTION definition, but with opposite signs. If needed, you can still add the combined REACTION incrementally, but I think the numerical simulation would be much more stable.

To remove a phase, for EQUILIBRIUM_PHASES, one way would be to set the target saturation index to be many orders of magnitude smaller with EQUILIBRIUIM_PHASES_MODIFY, it should not interfere with the numerical method.

Alternatively, if you are starting from scratch each time step, then you can simply define a new EQUILIBRIUM_PHASES with the phases that you want, where each phase has zero moles. Same with SOLID_SOLUTIONS.

If you switch to the differential REACTION approach, assuming the moles of the phase that you want to remove are zero, then you would define a new EQUILIBRIUM_PHASES with the phases you want to retain, but you would specify the number of moles of each from the last calculation for the cell. Same with SOLID_SOLUTIONS.

[sharmilan]

Dear Sir

Thank you for the reply, yes you are correct. The reaction in 100 steps at a particular time step is not possible in worker instance.

I have to run my mode for longer simulation times (100 years). as per the consideration of simulation time, I am able to set 12 hours as timestep for my simulation. Every 12 hours of the time step, I have to add nearly 1.5 mol of CO2 (R(n+1) - R(n)) to my chemical system and it makes trouble and convergence problem in PhrreqcRM when adding this larger mol amount. therefore, I tried to add this 1.5 mol of CO2 through REACTION keyword in 20 reaction steps with the setting of INCREMENTAL_REACTIONS true. In the InitialPhreeqc instance, this process is working well and include all this 1.5 mol at end of the reaction steps. but In worker, only one reaction step solute mol (1.5/20) added to the chemical system.

Therefore, is there any method to add 1.5 mol of CO2 in every time step into the system? otherwise, do I have to reduce my time step to 1/20 times?

Thank you
Sharmilan.S

[dlparkhurst]

At this point, I think you want to have some successful simulations, so you can see if your conceptual model is plausible. If using small time steps works, then run for a few days or months and see if problems persist. Otherwise, you can work on optimization after some initial testing.

Adding 1.5 moles of CO2 to a cell every 12 hours for 100 years does not make sense to me. I think that would be 100s of kg of carbon per 1 liter of water, but it is your simulation.

You could loop through RunCells 20 times for each transport step with REACTION defined as 1/20th of the total amount. However, the amount of time in the chemistry calculation would be the same as using a time step of 1/20th. Assuming the transport calculation is not prohibitive, you would be better solving transport (and chemistry) at the smaller time step, which should be more accurate, with little cost in execution time. RunCells would operate in parallel either way.

You could also run InitialPhreeqc-like calculations in the individual workers. So, you could do a multistep REACTION calculation. Using OpenMP, you can find the IPhreeqc instances and use IPhreeqc methods. Some example code is at https://phreeqcusers.org/index.php/topic,1686.msg5615.html#msg5615. However, this approach could be slower than running at a smaller time step. You would not have parallelization, so I think you would be worse off than just using a small time step and maximizing the parallel capabilities of PhreeqcRM.