RM_FindComponents and RM_GetComponent

[Yongqiang]

Hi experts,

I used RM_FindComponents to get the number of components in PHREEQCRM and RM_GetComponent to retrieve the names of components.
I found in my system that the number of components is 12. However, I retrieved only 11 names for components by RM_GetComponent. Could you please give me some advice about this discrepancy?

Thanks,
Yongqiang

[dlparkhurst]

I'm not sure what you did.

I'm assuming you are using Fortran. You should call RM_FindComponents after you have used RM_RunFile (or RM_RunString) for the InitialPhreeqc instance to define all of the reactants and solutions you are going to use in your calculation, and after you optionally set RM_SetComponentH2O.

Then you should use RM_GetComponents(id, components), where components is an array defined as follows:

Code: [Select]

character(len=:),   dimension(:), allocatable :: components

If you run RM_FindComponents, immediately followed by RM_GetComponents, I don't know how you would get different numbers of components. If you still get different results, I would need to see more of your input file and code.

[Yongqiang]

Hi David,

Please see the following code:
The code is in Matlab to call a compiled phreeqcrm library in the windows operating system.

Code: [Select]

clc
clear
nxyz = 40
loadlibrary('./bin/PhreeqcRMd','./include/RM_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')

ic = NaN(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

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

x6 = calllib('PhreeqcRMd','RM_SetComponentH2O',id, 0);

ncomps = calllib('PhreeqcRMd','RM_FindComponents',id);

% Parameters
bufferSize = 200;
% Preallocate a cell array to store components
components = cell(1, ncomps);
% Loop over components
for i = 1:ncomps
    disp(i);
    buffer = blanks(bufferSize);
    % Call the library function to get the component name
    [x6, buffer] = calllib('PhreeqcRMd', 'RM_GetComponent', id, i, buffer, bufferSize);
    disp(buffer)   
    % Store the modified buffer in the ith cell
    components{i} = strtrim(buffer); % Remove trailing spaces
    disp(i);
end


x998 = calllib('PhreeqcRMd','RM_CloseFiles',id)
x999 = calllib('PhreeqcRMd','RM_Destroy',id)
unloadlibrary PhreeqcRMd


[Yongqiang]

I updated the following code (C library is 0 based array, i corrected it in the following loop). A duplicate "charge" component is observed in the return values.

Code: [Select]

% Loop over components
for i = 1:ncomps
    disp(i);
    buffer = blanks(bufferSize);
    % Call the library function to get the component name
    [x6, buffer] = calllib('PhreeqcRMd', 'RM_GetComponent', id, i-1, buffer, bufferSize);
    disp(buffer)   
    % Store the modified buffer in the ith cell
    components{i} = strtrim(buffer); % Remove trailing spaces
    disp(i);
end

Code: [Select]

components

components =

  1?12 cell array

  Columns 1 through 11

    {'H'}    {'O'}    {'Charge'}    {'Al'}    {'C'}    {'Ca'}    {'Charge'}    {'Cl'}    {'Fe'}    {'Mg'}    {'Na'}

  Column 12

    {'Si'}

[dlparkhurst]

Let me see your input file.


Did you define an "element" named "Charge" in your input file?

[Yongqiang]

Hi David,

Please see the following input file:

Code: [Select]

KNOBS
-iterations 500
-convergence_tolerance 1e-8
-step_size 20
END

SOLUTION 0 formation water
     units   mol/kgw
     Na 0.5
     Cl 0.5
EQUILIBRIUM_PHASES 0
    CO2(g) 1.0 100
SAVE solution 0
END

SOLUTION 1 formation water
     units   mol/kgw
     #Na 0.00001
     #Cl 0.00001

EQUILIBRIUM_PHASES 1
    Albite 0 1499
    Anorthite 0 6347
    Forsterite 0 12086
    Fayalite 0 2746
    Diopside 0 295
SAVE SOLUTION 1
END

RATES

###########
#Fayalite
###########
Fayalite
# from Palandri and Kharaka 2004
# experimental condition range T=8-25C, pH=1.3-11

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=5.48E+11
12 E1=94400
13 n1=0
20 rem neutral solution parameters
21 a2=5.48E+02
22 E2=94400
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Fayalite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

############
#forsterite
############
Forsterite
# from Palandri and Kharaka 2004
# experimental condition range T=25-65C, pH=2-12

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=8.38E+04
12 E1=67206
13 n1=0.470
20 rem neutral solution parameters
21 a2=1.58E+03
22 E2=79000
30 rem base solution parameters
31 a3=1.00E-07
32 E3=56637
33 n2=-0.600
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("forsterite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

##############
#diopside
#############
Diopside
# from Palandri and Kharaka 2004
# experimental condition range T=8-90C, pH=1-6

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=3.00E+10
12 E1=96100
13 n1=0.710
20 rem neutral solution parameters
21 a2=1.00E-04
22 E2=40600
30 rem base solution parameters
31 a3=0
32 E3=0
33 n2=0
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("diopside")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

###########
#Albite
###########
Albite
-start
1 rem unit should be mol,Liter-1 and second-1
2 rem parm(1) is surface area in the unit of m2/L
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals M0 is the initial moles of minerals
5 rem parm(2) is a correction factor for fields rate relative to lab rate
10 rem acid solution parameters
11 a1=1.45E+0
12 E1=58400
13 n1=0.335
20 rem neutral solution parameters
21 a2=4.97E-10
22 E2=57000
30 rem base solution parameters
31 a3=7.15E-01
32 E3=55500
33 n2=0.317
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Albite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
110 rem do not dissolve more minerals than present
115 if (moles>M) then moles=M
200 save moles
-end

##############
#anorthite
##############
Anorthite
# from Palandri and Kharaka 2004
# experimental condition range T=25-95C, pH=2-10.2

-start
1 rem unit should be mol,kgw-1 and second-1
2 rem parm(1) is surface area in the unit of m2/kgw
3 rem calculation of surface area can be found in the note
4 rem M is current moles of minerals. M0 is the initial moles of minerals
5 rem parm(2) is a correction factor
10 rem acid solution parameters
11 a1=2.58E-01
12 E1=16601
13 n1=1.411
20 rem neutral solution parameters
21 a2=1.00E-06
22 E2=17821
30 rem base solution parameters
31 a3=1.00E-22
32 E3=18150
33 n2=-1.767
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("anorthite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("H+")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*parm(2)
100 moles= rate*Time
200 save moles
-end

KINETICS 0 define kinetic parameters
Fayalite
     -m 2.746#24.9
     -m0 2.746#24.9
     -parms 0.002  1


Forsterite
     -m 12.086
     -m0 12.086
     -parms 0.002  1
     
Albite
     -m 14.99
     -m0 14.99
     -parms 0.002  1


Anorthite
     -m 6.347
     -m0 6.347
     -parms 0.002  1
     
Diopside     
     -m 29.5
     -m0 29.5
     -parms 0.002  1   

END


SELECTED_OUTPUT
-reset false
-high_precision true
#-molalities Na+ Cl- H2O
USER_PUNCH
    -head rho pH molarity delta1 delta2 delta3 delta4 carbonates SI water
    10 PUNCH RHO
    20 PUNCH -LA("H+")
    30 PUNCH MOL("Fe+2")
    40 PUNCH MOL("Mg+2")
    50 PUNCH MOL("Ca+2")
    60 PUNCH KIN_DELTA("Fayalite")
    70 PUNCH KIN_DELTA("Forsterite")
    80 PUNCH KIN_DELTA("Diopside")
    90 PUNCH KIN_DELTA("Anorthite")
    100 PUNCH MOL("MgCO3")+MOL("FeCO3")+MOL("CaCO3")
    110 PUNCH SI("Calcite")
    120 PUNCH MOL("H2O")
END


[dlparkhurst]

That looks pretty normal.

My only other thought is that you did something like the following. If so, use SetComponentH2O only once, and before any call to FindComponents.

Code: [Select]

status = phreeqc_rm.SetComponentH2O(true);
ncomps = phreeqc_rm.FindComponents();
status = phreeqc_rm.SetComponentH2O(false);
ncomps = phreeqc_rm.FindComponents();


[Yongqiang]

Hi David,

Thanks for your insights!
You help me solve this problem. I did use duplicate findcomponents method. I removed the first one and everything is okay now.
I noticed a warning when running the code as the following. Do you have any idea what it is going on?

Thanks,

Code: [Select]

Warning: The function 'RM_BmiInitialize' was not found in the library
> In loadlibrary
In li (line 4)

[dlparkhurst]

Are you using Python or Fortran or C? Are you using a downloaded Python library? Or are you compiling a library?

My guess is that you the library you are using is not compiled with the USE_YAML definition. There may be a call to initialize with YAML that is not correctly #ifdef'd, which might cause that warning.

It would be best to compile with YAML, although getting YAML working can be an issue.

If you haven't used YAML and the library is working, I'd say don't worry about the warning message.

[Yongqiang]

Thanks, David.

Q:Are you using Python or Fortran or C? Are you using a downloaded Python library? Or are you compiling a library?
A:I am using the C methods in matlab. It is not a python anaconda installation. It is a direct compilation with Visual Studio on Windows operating system.

I found that the YAML method was introduced in the new version of PHREEQC. I am trying to figure it out what the benefits are to introduce YAML function. Is it for machine learning or managing large data?

Thanks,

[dlparkhurst]

No, the latest version of PhreeqcRM has a BMI (Basic Model Interface) that is supposed to simplify the process of coupling diverse models. In theory, each model uses the same methods. At the base level, there are only three methods, Initialize, Update (time step), and Finalize. Data are set with the SetValue method and retrieved with the GetValue method.

I did enough to satisfy the BMI standard, which I don't think really add that much value except for a couple features.

One is Initialize, which accepts a YAML file as the argument. There are a set of methods defined in the module YAMLPhreeqcRM that allow you to write initial conditions to a YAML file. So, you can define the number of cells, units, the initial temperature and pressure, solution and reactant numbers for model cells (InitialPhreeqc2Module), input file names, and other initial conditions and setting and write them to a YAML file. This file (along with other input and database files) would serve as an archive of your run. If you were to write a preprocessor for your model, you could use the YAML format to save the conditions defined in the preprocessor.

The second advantage of BMI is that there is a simpler way to retrieve variables with GetValue("varname", ...). A BMIPhreeqcRM method AddOutputVars lets you choose a set of variables that can be retrieved with GetValue. The method GetOutputVarNames provides a list of all variables that can be retrieved with GetValue. The bottom line is that you may not have to use SELECTED_OUTPUT/USER_PUNCH to extract the output variables that you need.

I have had trouble linking to some of the latest YAML releases. I use an older version (0.6.3) of YAML. But, if you succeed, everything you have done should be backward compatible. You can still use PhreeqcRM exactly as you have, except that you have access to the new BMI methods.

[Yongqiang]

Thanks for the awesome work, David.
This sounds fantastic. This new feature is worth learning to make the coupling simpler.

Cheers,