Obtain selected output from PhreeqcRM

[Ellis]

Hello everyone,

I am trying to obtain the selected_output from PhreeqcRM. I tried the following:

Code: [Select]

output = phreeqc_rm.GetSelectedOutput(true);
While the above outputs the species/components that I have specified in both the SELECTED_OUTPUT and USER_PUNCH data blocks, rather than an array of size m x n (where m = time_step and n is the number of specie/components) as I would have expected, what I get instead is an array of size k x n (where k = number of cells).

I would appreciate if someone can please kindly suggest to me what to do to get the correct output.

Thank you in anticipation.
Ellis

[dlparkhurst]

Although your syntax does not look correct, the results that you describe are correct. The C++ method will return a standard vector of size columns * number of cells (in sequence by row, that is the columns for row 1 followed by row 2, etc). The results will be for the current time step, logically following a call to RunCells. In your case, it sounds like this array would be columns of components and species and a row for each cell.

If you want results for all time steps and all cells, you will need to accumulate the results from each time step in other storage space, for example a map std::map<double, std::vector<double>, where the index of the map is the time and the std::vector is the result GetSelectedOutput. Alternatively, you can write the results to file to save the results for each time step.

"m x n (where m = time_step and n is the number of specie/components) as I would have expected". There is a set of species/components for each cell of the model, so are you expecting results from just one cell? Further, the method simply does not accumulate results over different time steps; it only returns values for the current time step.


The syntax for the call to GetSelectedOutput is as follows, where status is an int:

Code: [Select]


std::vector<double> so;
status = phreeqc_rm.GetSelectedOutput(so);
   

Size of the vector is set to col*nxyz, where col is the number of columns in the selected-output definition (GetSelectedOutputColumnCount), and nxyz is the number of grid cells in the user's model. (as described at https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm-html/class_phreeqc_r_m.html#a6179a2e85f8d5ac45174c097ac585dd0https://water.usgs.gov/water-resources/software/PHREEQC/documentation/phreeqcrm-html/class_phreeqc_r_m.html#a6179a2e85f8d5ac45174c097ac585dd0.)

A more complete example that considers the possibility of multiple selected output definitions is included in advection_cpp.cpp in the Tests directory of the distribution.

[Ellis]

Thank you for your response. I am actually using PhreeqcMatlab (see: https://github.com/simulkade/PhreeqcMatlab) as I do not understand C++.

"so are you expecting results from just one cell?". Yes, what I want is the selected_output results for the entire time steps of the last cell – just like I would normally get with Phreeqc Interactive/PhreeqcCOM.

The equivalent MATLAB syntax (courtesy of PhreecMatlab) for GetSelectedOutput is:

Code: [Select]

col = phreeqc_rm.RM_GetSelectedOutputColumnCount();
so  = zeros(col*nxyz, 1);
so  = reshape(so, nxyz, col);
[status, so] = phreeqc_rm.RM_GetSelectedOutput(so);
Thus, if you could give me the C++ syntax (or MATLAB, if possible) to obtain the accumulated results for the entire time steps of the last cell, then I should (or would) be able to find the equivalent MATALB code. Or as in your prescribed alternative, I would need the syntax to write the results to file (I would suppose this should be simple just like in COM?).

[dlparkhurst]

In C++ "so" is treated as a 1 dimensional array where the rows start at so[0], so[nxyz], etc. Looks like you are making it a two dimensional array, so it will be slightly different, and you may simplify not to use the temporary storage.

Code: [Select]

std::vector< std::vector<double> > all_times; # persistent storage
...
std::vector<double> so;
status = phreeqc_rm.GetSelectedOutput(so);
std::vector<double> last_cell;                          // temporary storage for last cell selected output
for (int i = 0; i < col; i++)
{
last_cell.push_back(so[(nxyz-1)*col + i]);  // Data for the last cell
}
all_times.pushback(last_cell);                            // transfer to persistent storage for all time planes

Alternatively, you can write so or all_times to file using Matlab methods.

[Ellis]

Thank you once again for your response. I have been able to extract the results for the last cell.

I then tried the code using the calcite dissolution example you described here: https://phreeqcusers.org/index.php/topic,1827.0.html except that I used my own transport solver, and I modified a few parameters (see attached).

The results for concentrations of Cl, Ca, (and C) agree (fairly) with what I get using Phreeqc interactive, whereas, there is a large discrepancy with what I get with moles of dissolved Calcite, and pH (see attached). Could my transport solver be the cause of this large discrepancy?

Also, you’ll notice that at step 0, M = 1 for Phreeqc interactive, but M = 0 in my code. I’ve tried to figure out why this is so but couldn’t. Same way I tried to figure out why I couldn’t get my ‘step’ to be computed incrementally. I really will appreciate your thoughts…

[dlparkhurst]

The "Calcite" columns is the value of M at each time, so PHREEQC has about 1.0 mole and your code has about 0.2 mole. I think the difference is a units conversion done by PhreeqcRM.

My guess is that you set a porosity of 0.2 and SetUnitsPPassemblage of 1, which would set the moles of calcite in the cell to be 1 (M in EQUILIBRIUM_PHASES definition) times porosity times representative volume (1 by default). So, the moles of calcite are different between the two simulations, as is the volume of water.  I think if if you use a representative volume of 5 L you would have 1 L of water and 1 mole of calcite and a porosity of 0.2. (Or you could use a porosity of 1 for PhreeqcRM, but not for your transport.)

[Ellis]

Thank you, Mr Parkhurst. I think I understand this better now.

However, I have been trying to play around with the local porosities of each cell. So, for instance, I assumed a particular cell is 100% solid, which implies porosity = 0 for that cell. When I run the code, I get an “unknown exception” error, which happens just before computing the initial boundary concentrations. I know this issue occurs because of the zero-porosity specified for one cell – as the code runs when I increase the porosity for that cell. So, please how can I overcome this? What should happen to the representative volume for this cell, and what kind of adjustment do I need to make the code to run?

[dlparkhurst]

You can use the saturation to determine whether a cell will be run or not. A zero saturation will cause PhreeqcRM to skip the calculation for that cell. Even so, I'm not sure you can set a zero porosity. You should probably use a small, nonzero number and set the saturation to zero.

[dlparkhurst]

If a transport-model cell is to have zero porosity throughout the calculation (an inactive cell), then you can use CreateMapping to avoid any chemistry calculations for that model cell. That transport-model cell will simply not have a chemistry cell equivalent, and its concentrations will be populated with missing values when you call GetConcentrations.

[Ellis]

Well, my intention is not to leave it that way all through the calculation. I intend to adjust/update porosity at every time step (if dissolution occurs). From what I have gathered from the forum, that can be made possible (for multicomponent diffusion) with CHANGE_POR. For now, I’m going one step at a time, so that I know where and how to fix an issue, if any.

Having said these, is there another way to proceed, bearing in mind that my zero porosity cell is not to be an “inactive cell”?

[Ellis]

Ok, I didn't see your penultimate response. My apologies. I think that setting the porosity to a small nonzero value might be a better approach. Let me try that out and see what happens.

PS: I already know about setting the saturation to zero if I do not want reaction calculations in a particular cell. Thanks to your help in the forum.