Initialization fails using PHREEQCRM with Matlab

[Yongqiang]

I used the following code in Matlab. x5 returns IRM_FAIL. Anyone can help with this issue?

Thanks,

Code: [Select]

clc
clear

nxyz = 40;
loadlibrary('./bin/PhreeqcRMd','./include/RM_interface_C.h');%RM_interface_C.hBMI_interface_C.h
id = calllib('PhreeqcRMd','RM_Create',nxyz,1);
libfunctionsview('PhreeqcRMd');
%libfunctions('PhreeqcRMd')
[x0,s] = calllib('PhreeqcRMd','RM_SetFilePrefix',id, 'Logfile');
x1 = calllib('PhreeqcRMd','RM_OpenFiles',id);
x2 = calllib('PhreeqcRMd','RM_LoadDatabase',id, 'LLNL.DAT');
%tt = calllib('PhreeqcRMd','RM_FindComponents',id);
[x3,s] = calllib('PhreeqcRMd','RM_RunString',id,1,1,1,'END');
[x4,s] = calllib('PhreeqcRMd','RM_RunFile',id,1,1,1,'./solution.pqi');
x4 = calllib('PhreeqcRMd','RM_RunString',id, 1, 0, 1, 'DELETE; -all');

ic = ones(1,7*nxyz);
for i = 1:nxyz
    ic(i) = i;
    ic(i+nxyz) = i;
    ic(i+2*nxyz) = -1;
    ic(i+3*nxyz) = -1;
    ic(i+4*nxyz) = -1;
    ic(i+5*nxyz) = -1;
    ic(i+6*nxyz) = -1;
end
t1 = libpointer('int32Ptr', []);
t2 = libpointer('doublePtr', []);
[x5,s1,s2,s3] = calllib('PhreeqcRMd','RM_InitialPhreeqc2Module',id, ic, t1, t2);

[dlparkhurst]

I don't use Matlab regularly, so I don't really want to debug a script. Here is a previous post with a script that should run. Perhaps that will get you started.

https://phreeqcusers.org/index.php/topic,2564.msg9640.html#msg9640

[Yongqiang]

Thanks, David.
The library loading and phreeqc instance creation are successful. Only the following code fails, could you please help me with it? Thank you!

Code: [Select]

[x5,s1,s2,s3] = calllib('PhreeqcRMd','RM_InitialPhreeqc2Module',id, ic, t1, t2);

[dlparkhurst]

Set t1 to -1 and t2 to 1.0, same dimension as ic1.

[Yongqiang]

Thank you so much, David.
The following is my full code. However, I tried the comments in above. It forced quit the Matlab program.

Code: [Select]

clc
clear

nxyz = 40;
loadlibrary('./bin/PhreeqcRMd','./include/RM_interface_C.h');%RM_interface_C.hBMI_interface_C.h
id = calllib('PhreeqcRMd','RM_Create',nxyz,1);

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

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

% Set PhreeqcRM time units
status = calllib('PhreeqcRMd', 'RM_SetTimeConversion', id, 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(nxyz, 1, 'double');
status = calllib('PhreeqcRMd', 'RM_SetRepresentativeVolume', id, 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', id, por);
% Set PhreeqcRM saturation
sat = ones(nxyz, 1, 'double');
status1 = calllib('PhreeqcRMd', 'RM_SetSaturation', id, sat);
% Disable the detailed printing of chemical information
status2 = calllib('PhreeqcRMd', 'RM_SetPrintChemistryOn', id, 0, 1, 0);
%% Get cell numbers
nchem = calllib('PhreeqcRMd', 'RM_GetChemistryCellCount', id);

libfunctionsview('PhreeqcRMd');
%libfunctions('PhreeqcRMd')
[x0,s] = calllib('PhreeqcRMd','RM_SetFilePrefix',id, 'Logfile');
x1 = calllib('PhreeqcRMd','RM_OpenFiles',id);
x2 = calllib('PhreeqcRMd','RM_LoadDatabase',id, 'LLNL.DAT');
%tt = calllib('PhreeqcRMd','RM_FindComponents',id);
[x3,s] = calllib('PhreeqcRMd','RM_RunString',id,1,1,1,'END');
[x4,s] = calllib('PhreeqcRMd','RM_RunFile',id,1,1,1,'./solution.pqi');
x4 = calllib('PhreeqcRMd','RM_RunString',id, 1, 0, 1, 'DELETE; -all');
%%
ic = ones(7*nxyz,1);
ic1 = ones(7*nxyz,1);
f = ones(7*nxyz,1);
for i = 1:nxyz
    ic(i) = i;
    ic(i+nxyz) = i;
    ic(i+2*nxyz) = -1;
    ic(i+3*nxyz) = -1;
    ic(i+4*nxyz) = -1;
    ic(i+5*nxyz) = -1;
    ic(i+6*nxyz) = -1;
   
    ic1(i) = -1;
    ic1(i+nxyz) = -1;
    ic1(i+2*nxyz) = -1;
    ic1(i+3*nxyz) = -1;
    ic1(i+4*nxyz) = -1;
    ic1(i+5*nxyz) = -1;
    ic1(i+6*nxyz) = -1;

    f(i) = 1.0;
    f(i+nxyz) = 1.0;
    f(i+2*nxyz) = 1.0;
    f(i+3*nxyz) = 1.0;
    f(i+4*nxyz) = 1.0;
    f(i+5*nxyz) = 1.0;
    f(i+6*nxyz) = 1.0;
end
x6 = calllib('PhreeqcRMd','RM_SetComponentH2O',id, 0);
%%
%t1 = libpointer('int32Ptr', []);
%t2 = libpointer('doublePtr', []);
[x5,s1,s2,s3] = calllib('PhreeqcRMd','RM_InitialPhreeqc2Module',id, ic, ic1, f);%
I tried to run section by section and identified the problem is still the RM_InitialPhreeqc2Module command. The last line:

Code: [Select]

[x5,s1,s2,s3] = calllib('PhreeqcRMd','RM_InitialPhreeqc2Module',id, ic, ic1, f);
The ic1 and f were set with the following lines:

Code: [Select]

ic = ones(7*nxyz,1);
ic1 = ones(7*nxyz,1);
f = ones(7*nxyz,1);
for i = 1:nxyz
    ic(i) = i;
    ic(i+nxyz) = i;
    ic(i+2*nxyz) = -1;
    ic(i+3*nxyz) = -1;
    ic(i+4*nxyz) = -1;
    ic(i+5*nxyz) = -1;
    ic(i+6*nxyz) = -1;
   
    ic1(i) = -1;
    ic1(i+nxyz) = -1;
    ic1(i+2*nxyz) = -1;
    ic1(i+3*nxyz) = -1;
    ic1(i+4*nxyz) = -1;
    ic1(i+5*nxyz) = -1;
    ic1(i+6*nxyz) = -1;

    f(i) = 1.0;
    f(i+nxyz) = 1.0;
    f(i+2*nxyz) = 1.0;
    f(i+3*nxyz) = 1.0;
    f(i+4*nxyz) = 1.0;
    f(i+5*nxyz) = 1.0;
    f(i+6*nxyz) = 1.0;
endIs there any way to fix it?

Thanks!

[dlparkhurst]

Is this right? Do you have 40 different solutions and equilibrium_phases defined in the input file? Sure you don't mean 1 instead of i?

Code: [Select]

    ic(i) = i;
    ic(i+nxyz) = i;


[Yongqiang]

Ah, that makes sense. How stupid the mistake is! Thank you so much, David!!!!!