PhreeqcUsers Discussion Forum
Beginners => PHREEQC basics => Topic started by: Minjeong on March 16, 2020, 08:55:05 AM
-
Good day!
I'm currently running Iphreeqc with MATLAB.
I compared the results of two code blocks and they showed different results for charge balance though I expected they showed the same results.
The difference between them is initialization of the Iphreeqc object before the command SOLUTION_MODIFY.
Chemical.DB='phreeqc.dat';
Chemical.CH={'Ca','Cl','K','N','Na','O(0)'}';
PStr1=sprintf(['SOLUTION 1\n', ...
'-units mol/L\n', ...
'pH 7 charge\n', ...
'Ca 0.0006\n', ...
'Cl 0.0012\n', ...
'O(0) 1e-5\n']);
PStr2=sprintf(['SOLUTION_MODIFY 1\n',...
'-totals\n',...
'Ca 0.005\n',...
'Cl 0.003\n',...
'K 0.02\n',...
'N 0\n',...
'Na 0\n',...
'O(0) 1.000066756e-05\n',...
'temp 32.400000\n']);
STR = sprintf(['SELECTED_OUTPUT\n', ...
'-reset false\n', ...
'-time\n',...
'-high_precision true\n', ...
'-solution\n',...
'-water\n',...
'-pH\n',...
'-pe\n',...
'-charge_balance\n',...
'-percent_error\n',...
'-equilibrium_phases Calcite\n',...
'-saturation_indices Calcite\n',...
'-molalities NaX CaX X-\n',...
'USER_PUNCH\n', ...
['headings', ' density', cell2mat(strcat({' '},Chemical.CH')), '\n'], ...
'-start\n', ...
['10 punch', ' rho,', cell2mat(strcat({' totmol("'},Chemical.CH',{'"),'})), '\n'], ...
'-end\n', ...
'END\n']);
[b]%% case1) using a single Iphreeqc obeject w/o INCLUDE$ dumpfile[/b]
iphreeqc = actxserver('IPhreeqcCOM.Object');
iphreeqc.LoadDatabase(fullfile(cd, 'Phreeqc_DB', Chemical.DB));
iphreeqc.DumpFileOn = true;
iphreeqc.RunString(PStr1); %SOLUTION
iphreeqc.RunString(STR); %SELECTED_OUT
iphreeqc.AccumulateLine('RUN_CELLS');
iphreeqc.AccumulateLine('-cells 1');
iphreeqc.AccumulateLine('-time_step 0');
iphreeqc.AccumulateLine('DUMP')
iphreeqc.AccumulateLine('-file dump_initial')
iphreeqc.AccumulateLine('-all')
iphreeqc.RunAccumulated;
output_inital=iphreeqc.GetSelectedOutputArray;
iphreeqc.ClearAccumulatedLines;
iphreeqc.AccumulateLine(PStr2); %SOLUTION_MODIFY
iphreeqc.AccumulateLine('RUN_CELLS');
iphreeqc.AccumulateLine('-cells 1');
iphreeqc.AccumulateLine('DUMP')
iphreeqc.AccumulateLine('-file dump_output1')
iphreeqc.AccumulateLine('-all')
iphreeqc.RunAccumulated;
output1=iphreeqc.GetSelectedOutputArray;
[b]%% case2) reset iphreeqc object w/ INCLUDE$ dumpfile[/b]
iphreeqc.delete;
iphreeqc = actxserver('IPhreeqcCOM.Object');
iphreeqc.LoadDatabase(fullfile(cd, 'Phreeqc_DB', Chemical.DB));
iphreeqc.DumpFileOn = true;
iphreeqc.RunString('INCLUDE$ dump_initial') %INITIALIZE the solution compositions
iphreeqc.AccumulateLine(PStr2); %SOLUTION_MODIFY
iphreeqc.AccumulateLine(STR); %SELECTED_OUT
iphreeqc.AccumulateLine('RUN_CELLS');
iphreeqc.AccumulateLine('-cells 1');
iphreeqc.AccumulateLine('DUMP')
iphreeqc.AccumulateLine('-file dump_output2')
iphreeqc.AccumulateLine('-all')
iphreeqc.RunAccumulated;
output2=iphreeqc.GetSelectedOutputArray;
The results of two cases (dump_output1 and dump_output2) are different for charge balance (plus, some small digits of H(0), pe, etc).
I would like to know how to make them to show the totally same results.
Thanks in advance.
Sincerely,
Minjeong
-
I think you are complaining about the difference in charge balance between +1e-17 and -1e-17. Sorry, but you will not get better resolution than that. In fact, the default convergence tolerance is 1e-8 times the total concentration for each element and 1e-8 for charge balance, so the result could be worse. You can decrease the convergence tolerance in KNOBS, say to 1e-10 or 1e-12, to ensure a bit more precision. However, an estimate of roundoff error is about 1 in the 16th decimal digit for a typical 64-bit computer; so 1e-16 is the minimum possible value for convergence tolerance, but realistically you are likely to have convergence problems below 1e-12 for convergence tolerance.