Using selected output in phreeqcRM

[Wonwoo Yoon]

Hello all,

I'm currently working on coupling phreeqcRM and COMSOL.

I'm okay with other parts of the code, and I successfully benchmarked Example 11.

However, I'm still having a hard time utilizing SELECTED_OUTPUT in phreeqcRM and want to get some help with this.

Since I'm using MATLAB for the data transfer between phreeqcRM and COMSOL, I load the C++ library of phreeqcRM in MATLAB.

In the code below,

status = calllib('phreeqcRMd', 'function_name', input_arguement_1, input_arguement_2);

is equivalent to

status = function_name(input_argument_1, input_argument_2).

The following is the initial part of my code.

Code: [Select]

% PhreeqcRM node generation
phr_model = calllib('PhreeqcRMd', 'RM_Create', node, -1);

% Set PhreeqcRM properties
status = calllib('PhreeqcRMd', 'RM_SetErrorOn', phr_model, 1); %#ok<*NASGU>
status = calllib('PhreeqcRMd', 'RM_SetErrorHandlerMode', phr_model, 2);
status = calllib('PhreeqcRMd', 'RM_SetComponentH2O', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetRebalanceFraction', phr_model, 0.5);
status = calllib('PhreeqcRMd', 'RM_UseSolutionDensityVolume', phr_model, 0);
status = calllib('PhreeqcRMd', 'RM_SetPartitionUZSolids', phr_model, 0);
status = calllib('PhreeqcRMd', 'RM_SetFilePrefix', phr_model, 'react_c');

% Set PhreeqcRM units
status = calllib('PhreeqcRMd', 'RM_SetUnitsSolution', phr_model, 2);
status = calllib('PhreeqcRMd', 'RM_SetUnitsPPassemblage', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetUnitsExchange', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetUnitsSurface', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetUnitsGasPhase', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetUnitsSSassemblage', phr_model, 1);
status = calllib('PhreeqcRMd', 'RM_SetUnitsKinetics', phr_model, 1);

% Set PhreeqcRM time units
status = calllib('PhreeqcRMd', 'RM_SetTimeConversion', phr_model, 1); % number x represents 1 second = x timestep

% Set PhreeqcRM representative volume, set a large rv when the pore volumn * saturation is low
rv = ones(node, 1, 'double');
status = calllib('PhreeqcRMd', 'RM_SetRepresentativeVolume', phr_model, rv);

% Set PhreeqcRM porosity
porosity_index = find(strcmp(transport_data.expr, 'dl.epsilon'));
por = transport_data.(['d',num2str(porosity_index)])';

status = calllib('PhreeqcRMd', 'RM_SetPorosity', phr_model, por);

% Set PhreeqcRM saturation
sat = ones(node, 1, 'double');
status = calllib('PhreeqcRMd', 'RM_SetSaturation', phr_model, sat);

% Disable the detailed printing of chemical information
status = calllib('PhreeqcRMd', 'RM_SetPrintChemistryOn', phr_model, 0, 1, 0);

% Get cell numbers
nchem = calllib('PhreeqcRMd', 'RM_GetChemistryCellCount', phr_model);

% Load database
status = calllib('PhreeqcRMd', 'RM_LoadDatabase', phr_model, 'C:\Program Files\USGS\IPhreeqcCOM 3.7.3-15968\database\phreeqc.dat');

% Run initial run for each chemical zone
% The initial condition for each chemical zone should be defined in 'initial.phr', before the model run.
status = calllib('PhreeqcRMd', 'RM_RunFile', phr_model, 1, 1, 1, 'initial.phr');
status = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, 'DELETE\n -all\n');

and the 'initial.phr' for my code is as below.

Code: [Select]

SOLUTION 1
units            mmol/kgw
temp             25.0
pH               7.0     charge
pe               12.5    O2(g)   -0.68
Na               1.0
K                0.2
N(5)             1.2
END

EXCHANGE 1
-equilibrate 1
X 0.0011
END

SOLUTION 2
units            mmol/kgw
temp             25.0
pH               7.0     charge
pe               12.5    O2(g)   -0.68
Ca               0.6
Cl               1.2
END

EXCHANGE 2
-equilibrate 2
X 0.0011
END

SELECTED_OUTPUT
-reset FALSE
-high_precision TRUE
USER_PUNCH
-headings CaX2 KX NaX
10 PUNCH MOL('CaX2') MOL('KX') MOL('NaX')
END

Actually, the code above worked really well without the SELECTED_OUTPUT part in 'initial.phr', but it always failed to run after I added the SELECTED_OUTPUT at line 'status = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, 'DELETE\n -all\n');'.

The error message in MATLAB was

Code: [Select]

Debug Assertion Failed!
C:\my_directory\PhreeqcRMd.dll
C:\phreeqcrm_source\phreeqcrm-3.7.3-15968-hotfix\src\IPhreeqcPhast\IPreeqc\IPhreeqc.cpp
Line: 369

Expression: (null)

Line 367-371 in IPhreeqc.cpp was as below:

Code: [Select]

int IPhreeqc::GetSelectedOutputCount(void)const
{
ASSERT(this->PhreeqcPtr->SelectedOutput_map.size() == this->SelectedOutputMap.size());
return (int) this->PhreeqcPtr->SelectedOutput_map.size();
}

I would greatly appreciate it if you provide a solution to this problem.

Thanks for reading.

Sincerely, Woowoo Yoon

[Wonwoo Yoon]

One mistake from my original text.

I imported C library of phreeqcRM to MATLAB not C++.

[dlparkhurst]

I don't see a problem with the code that you show so far, although I may have missed something.

You do need to distribute initial conditions (InitialPhreeqc2Module) and then run the cells (RunCells) before selected output is available to be retrieved. If you have have example 11 working, then try to follow the same sequence of calls.

[Wonwoo Yoon]

I found out that deleting the code

'status = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, 'DELETE\n -all\n');'

makes it work. However, as I understand, this code was just for clearing out (maybe initialization?) the module and utility before InitialPhreeqc2Module.

I'm not pretty sure why this is happening, but are there any further problems that I might encounter if I don't initialize the module and utility and just copy the conditions in Initialphreeqc to the module by using InitialPhreeqc2Module?

[dlparkhurst]

I don't see why that call should have caused a problem. I'll see if I can reproduce it. It is an "assert", which will only happen in the debug version; its a guess, but I think the release version would probably run OK.

The reason I do the DELETE is that when you run the initial.phr, it does create SOLUTIONS, EQUILIBRIUM_PHASES, etc in each of the IPhreeqc instances. So, if you define EQUILIBRIUM_PHASES 2 in initial.phr, for example, cell 2 will have an equilibrium phases definition (in whichever IPhreeqc instance that is assigned cell 2). Now, if you do InitialPhreeqc2Module, if you do specify an equilibrium phases for cell 2, it will overwrite the initial definition, and things will work correctly. The problem would be if you did not intend to have equilibrium phases for cell 2, the original definition would still exist, and cell 2 would in fact have an equilibrium phases definition.

[dlparkhurst]

Another option is to define SELECTED_OUTPUT and any thermodynamic database changes (SOLUTION_MASTER_SPECIES, SOLUTION_SPECIES, PHASES, etc) in a file that you run in all the workers, initial phreeqc, and utility instances, and then run any any SOLUTIONs, EQUILIBRIUM_PHASES, etc only in the InitialPhreeqc instance. Afterward you would use InititialPhreeqc2Module to populate the cells and there would be no danger of stray definitions.

But that requires two files and careful separation of the data blocks into the appropriate file. It seemed easier to clear out the definitions from the workers and utility.

[dlparkhurst]

I still don't understand why you had a problem. I would like to know if it runs if you choose 1 thread in the constructor RM_Create. Perhaps you had more threads than cells? But that should result in a different error message.

Also, reading your post more carefully, we plan to discontinue support for C. We were planning to support C++, Fortran, and Python.

[Wonwoo Yoon]

I think maybe the version of phreeqcRM is the problem.

I experienced the LNK2019 ERROR during the installation of phreeqcRM (openMP, .dll), which had been already discussed in this forum (https://phreeqcusers.org/index.php/topic,2166.msg7960.html#msg7960).

I downloaded the hotfix version from the above post, and maybe it could be a reason for this error.

I will try it with the .lib version of 3.7.3-15968 (uploaded in USGS), since the LNK2019 ERROR is known to occur only in Windows .dll version.

Additionally, I will have to figure out how to use the C++, Fortran, or Python version of phreeqcRM through MATLAB in that case. Thanks for your announcement.

[dlparkhurst]

I think we should continue off the forum. I will send you a message and we can continue the discussion.

[dlparkhurst]

The problem appears to be in the statement

Code: [Select]

status = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, 'DELETE\n -all\n');

Matlab does not interpret \n as the newline character. It interprets it as two characters: a backslash ('\') and an n ('n').

Here are a couple of alternatives. You can use a semicolon, which is interpreted as the end of a line by Phreeqc:

Code: [Select]

tatus = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, 'DELETE; -all');

Or you can use the Matlab newline function as follows:

Code: [Select]

status = calllib('PhreeqcRMd', 'RM_RunString', phr_model, 1, 0, 1, ['DELETE' newline '-all' newline]);


[dlparkhurst]

You should consider using PhreeqcMatlab as an easier interface from Matlab to PhreeqcRM.

https://www.mathworks.com/matlabcentral/fileexchange/99394-phreeqcmatlab