PhreeqcRM failed in the middle of the numerical simulation

[Zhaoyang]

Dear all,

I encountered a problem when running the model that couples PhreeqcRM with a transport model. In the middle of the numerical simulation (30350 out of 72000 steps), PhreeqcRM failed because of

"ERROR: RunString: bad allocation; Stopping; ERROR: Run Cells failed in worker 0 from an unhandled exception; ERROR: RunString: bad allocation; ERROR: PhreeqcRM failed; PhreeqcRM:: RunCells."

I checked the latest result from the coupling model before the error occurred. The result is quite normal before the model reaches the steady state. To fix this problem, I changed parameters for the transport model, including setting a finer mesh, increasing the dispersivity which makes the model more stable. However, all of things I have tried don't work. Does anyone have experience in solving this problem? Thanks for your help in advance.

Best regards,
Zhaoyang

[dlparkhurst]

Maybe you ran out of memory?  Otherwise, I really don't know.

[Zhaoyang]

Thanks for your reply.

To work it out, I have tried different ways including a finer mesh, a larger dispersivity coefficient and a smaller time step. These usually make the numerical model more stable, but the problem is still there. In addition, I notice that the simulation with the same parameters stopped at different steps. Therefore, I think you are right. I was wondering if I can increase the memory to solve this problem. It should be noted that the latest input from the model is reasonable even though the simulation failed in the middle.

Furthermore, I have another question about implementing warm start in my coupling model. I have 17141 nodes (61 in vertical direction and 281 in horizontal direction). I want to use warm start to save time. However, whatever concentrations I set to the transport model, the initial concentrations are always zero for all secies. I guess there are some problems in:

Code: [Select]

  ! Set printing of chemistry file to false
  status = RM_SetPrintChemistryOn(id, 0, 1, 0)  ! workers, initial_phreeqc, utility
  ! Load database
  status = RM_LoadDatabase(id, "phreeqc.dat")
  status = RM_RunFile(id, 1, 1, 1, "advect.pqi")
  ! Clear contents of workers and utility
  string = "DELETE; -all"
  status = RM_RunString(id, 1, 0, 1, string)  ! workers, initial_phreeqc, utility

Could you please give me some suggestions and explain a bit more about workers, initial_phreeqc and utility? Thanks for your help in advance.

By the way, I found that I can't upload files to this forum.

[dlparkhurst]

The error message that you included in your first post looks like you ran out of memory and it failed in an effort to allocate more memory. If you are running on Windows, you can use the TaskManager to look at the memory consumption of your program; on Linux you can use "top" command among others.

If it is memory, it could be your side or PhreeqcRM. PhreeqcRM does not have any memory leaks that we know of, but there are probably some ways you could continuously use more memory. Multiple calls to SaveState without calling StateDelete would use more and more memory. Running with a coarser grid will require less memory, so you could test whether you get farther in your simulation as you use a coarser grid (or run less model time with a finer grid). Using a computer with more memory is obviously a possibility.

I don't think that deleting all of the cells in the workers for a warm start is possible or effective. You will simply have to recreate the definitions for each cell, which should use the same amount of memory. Plus, you will be deleting the state of all of the reactants (like EQUILIBRIUM_PHASES) which you would have to somehow retrieve from a saved state.

You could probably delete all of the definitions in initial_phreeqc once you populated all of the cells in the worker instances (except for changing boundary condition information). The intention is that there are just a few definitions in the initial_phreeqc instance; an end-member case is that initial conditions are the same throughout your domain, in which case only one set of initial solution and reaction definitions in initial_phreeqc would suffice to populate all of the cells. In this case, there is not much memory associated with the initial_phreeqc instance. Depending on how you do it, you could possibly have separate definitions for each cell in the initial_phreeqc instance, in which case there could be a significant savings by deleting everything in initial_phreeqc. Usually, the utility instance is not used, so there usually is not a lot of memory associated with it, but, depending on the logic of your program, you might be able to delete everything in the utility instance.

Usually, I think the amount of memory required would be similar from one time step to the next in PhreeqcRM. So, I am not sure what would be using more and more memory as the program progresses.

[Zhaoyang]

By monitoring memory, I am sure the stop is because of running out of memory. Memory used by the simulation first increases and then decreases. After that, the memory used by the simulation increases again until the simulation fails.

To figure out this problem, I run the simulation without reactions (turning of PhreeqcRM). In this way, the simulation ends normally. In addition, memory used by the simulation is stable during the whole simulation. Do you have any suggestions? Thanks for your help in advance.

Best regards,
Zhaoyang

[dlparkhurst]

Limit the number of chemical components to those that you are really interested in (PHAST transports all the elements in SOLUTIONs, EQUILIBRIUM_PHASES, KINETICS, etc, although not the surface like Hfo or exchangers like X).

Use a coarser grid.

Use a computer with more memory.

[Ivan K.]

I'm trying to solve a reactive transport problem with 10500 cells. Due to the necessity of a sharp dissolution front tracking a small time step is used. Thus I have to make tens of thousands of time steps. In the Task manager I observe continuous growth of memory consumption of the program. Finally the calculation crashes due to running out of memory. 64 Gb of RAM are installed on my desktop. The problem vanishes if I turn off chemistry by just commenting one line of code with a call to RunCells(). In the latter case the memory consumption is stable along the simulation.
My PHREEQCRM version is 3.7.0-15749. I will try a newer one.

[Ivan K.]

I've just compiled and tried the last version PHREEQCRM3.7.3-15968 but the problem remains.

[Ivan K.]

It seems that some data is continuously stored by PHREEQCRM on each time step. The occupied memory is freed after the computation is finished and PHREEQCRM object is deleted. So, it probably is not a memory leak, but some excessive data is stored during the simulation.

[dlparkhurst]

I'm not sure what is using all the memory. We check for memory pretty carefully. I don't think 11000 cells should be a problem

Two unlikely possibilities are a memory checker that does not reuse freed memory or if you save many states (StateSave).

You can try eliminating all output -- selected output, chemistry print and dump, to see if a lot less memory is used.

You can also eliminate reaction--equilibrium phases, kinetics etc -- to see if they are using a lot of memory.

[Ivan K.]

David, thanks for the advice.
I don't save states at all.
But I've found the source of the problem - it is the SELECTED OUTPUT. I tried to turn it off completely and, yes, the memory consumption became stable during the simulation.
It seems that output is accumulated somewhere within PHREEQCRM. Is there any possibility to clean the corresponding data?
On each time step I obtain only a fixed-size vector<double> of selected output, thus I expected only this vector to be stored. I checked that the size of this vector is constant.
Actually I need Selected Output as I retrieve some data from the simulation using it, thus turning off output doesn't seem to be a good solution.

[dlparkhurst]

Are you getting successively greater memory usage at each iteration? I don't think that should happen, and if it does, it would be a bug. Normally, you would define SELECTED_OUTPUT/USER_PUNCH for the workers before you begin time stepping, and then that definition would continue through the simulation.

If you need the selected output, then I think you will just have to spend the memory to get it, but it should not increase unless you are adding additional SELECTED_OUTPUT definitions.

Is your selected output definition for a single time step really large?

[Ivan K.]

Yes, I'm getting successively greater memory usage at each time step. I see it continuously growing. SELECTED_OUTPUT is defined only once in the beginning before the time stepping starts. Every time step the used memory increases by 2-3 Mb. I managed to overcome the problem by demanding this output less frequently, not on every time step.
I'm using    phreeqc_rm->SetSelectedOutputOn(***) function to turn output on and off during the simulation.
Thus the program fits in memory. But it seems to be a bug.
The memory consumption for sure doesn't increase if I turn off selected output.

The demanded output is not large, I request pH,pe and one equlibrium phase.

[dlparkhurst]

What operating system are you using?

What compiler are you using?

What compiler options or Visual Studio configuration are you using?

[dlparkhurst]

Are you using MPI or OpenMP?

[Ivan K.]

I'm using Windows 10 OS. PHREEQCRM is compiled by Visual Studio 2017 compiler. I compile using the instructions given in README.TXT.
OpenMP is used.

[Ivan K.]

I made an important observation. The memory consumption strongly depends on the number of openMP threads. I ran my simulation for 50 time steps with different number of threads and got the following memory usage at the end:
#threads   RAM used, Mb
1              287
2              302
4              359
8              470
16            750
These values are approximate (I just observed the maximum values in the Task Manager), but the dependence is clear.

I also understood that for me less frequent output is not a good solution because I need at least to get equilibrium phases on every time step.

[dlparkhurst]

Thanks for the information. I don't think it is really a memory leak; we check that every night. My hypothesis is that it is something in the way the standard template classes are handled where memory is not cleaned up until later, but it could be something we programmed. We are starting to investigate.

[Ivan K.]

I also don't think that it's a memory leak which can be caught by valgrind. This memory is freed when the simulation is over but the software is still running. Thus, it is most likely freed when PHREEQCRM instance is deleted. So, this memory seems to be consumed by some internal data storage.

[dlparkhurst]

Scott tracked down the bug. A selected-output line was continually added to an array without ever clearing the the array. We'll try to get you a fix soon.

[Ivan K.]

This is good news! I'm looking forward to obtain the corrected version.

[charlton]

Hopefully this fixes the problem:

https://github.com/usgs-coupled/phreeqcrm/releases/tag/v3.7.6-rc1

Scott

[Zhaoyang]

Dear Scott,

Thank you so much for the update. This problem has annoyed me for a long time.

Best regards,
Zhaoyang

[Ivan K.]

Hopefully this fixes the problem:

https://github.com/usgs-coupled/phreeqcrm/releases/tag/v3.7.6-rc1

Scott

Dear Scott, I'm testing this version. It appears that the function SetSaturation is now called SetSaturationUser. Is that done intentionally?

[dlparkhurst]

The version Scott posted is actually a release candidate for the next release. There are many additions to PhreeqcRM, but in terms of usage, you have pointed out the main non-backward compatible change. I thought it was confusing to have SetSaturation and GetSaturation actually refer to different arrays. Same with SetDensity and GetDensity. Here is the new documentation:

-----------------
May 19, 2023
-----------------
PhreeqcRM:
Renamed GetDensity and related functions to GetDensityCalculated.
Renamed SetDensity and related functions to SetDensityUser.

Density is used to convert user-model concentrations to module solution definitions only if
the units of the user-model concentrations are specified to be parts per million. The
density specified by SetDensityUser is used by SetConcentrations to convert from per kg of
solution to per L of solution. For GetConcentrations, two options are available to convert
from module solutions to user-model concentrations, depending on the value used for the
method SetUseSolutionDensityVolume: (1) the module-calculated density is used to convert
from the calculated volume of solution to the mass (kg) of solution, or (2) the
user-specified value of density is used to make the conversion. Again, density is only used
if the user-model concentration units are ppm.

The change in method names is intended to emphasize the difference between the
user-specified densities and the module-calculated densities.

Renamed GetSaturation and related functions to GetSaturationCalculated.
Renamed SetSaturation and related functions to SetSaturationUser.

The values specified by SetSaturation are used to convert user-model concentrations to
module solution definitions. For SetConcentrations, the volume of solution is calculated to
be the user-specified saturation * porosity * representative volume. For GetConcentrations,
two options are available to determine the solution volume, depending on the value
specified for SetUseSolutionDensityVolume: (1) the solution volume is calculated by the
reaction module and used to convert to user-model concentrations, or (2) the solution
volume is again calculated by user-specified saturation * porosity * representative volume,
and those values are used to convert to user-model concentrations. In either case, the
values returned by GetSaturationCalculated are the calculated solution volume divided by
(porosity * representative volume).
   
The change in method names is intended to emphasize the difference between the
user-specified saturations and and the module-calculated saturations.