mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-16 00:28:23 +01:00
b8745514b6
2243d25 Merge commit '013c822f76e5dc2e4fc19e87c6e5777aea6151d2' c1af6f3 added newlines for CRAN 013c822 added newlines for CRAN e4bd9ba [phreeqc3] fixes -Wclass-memaccess warnings for CRAN 29f06d2 fixed alignment in Description of solution 09a2680 guarded write_banner with NO_UTF8_ENCODING 082edbb changed src/print.cpp back to windows-1252 encoding; updated check_utf.sh 8d7c1fc adding mcd_Jtot and mcd_Jconc 9f0f622 Merge branch 'master' of github.com:usgs-coupled/phreeqc3 1040066 Merge remote-tracking branch 'usgs-coupled/master' 2a94644 cleaned up to eliminate some prints 07a864d all jacobians are consistent. Looks pretty good. 56975a7 Saved surface for numerical derivatives df0d68b Runs all the test cases. Numerical derivatives work, but still some changes in residuals before and after jacobian calculations. 6bd936e Fixed numerical derivative (non-pitzer) 0dde2b0 removed comments aef51fa Finally have derivatives right, I think 20281a0 always reset gases 13ec2fc best I could do for H2S while maintaining old tests. Used INCREMENTAL reactions 8be1ba8 revised jacobian_pz with new logic. Works with fixed_pressure examples H2S, H2S_pz, H2S_pz_appt, H2S_NaCl_Na2SO4. 71cf2a9 still produces different residuals 9022ded Tony H2S. Amm.dat, phreeqc.dat, pitzer.dat, utf8, updated test cases cb1f9af Finished up C, Fortran, documentation. Need to check DOxygen 9dad447 Merge remote-tracking branch 'origin/master' into state d647eec Added StateSave, StateApply, StateDelete with documentation for C++. Need testing, Fortran, and C 48cb5e8 Including OH- in converting units. Revised calculated density for H+ and OH-. Makes a difference in several test cases. Removed timing at end of .out in test cases. Checking in all test cases and selected output. 47e1ce5 added OH in density iteration calculation, test case NaOH_density 4aefb06 allow Fe(+3), equivalent to Fe(3), in TOT and TOTMOL. Previously fixed in SELECTED_OUTPUT -total bea0ad1 unused variable, punch Fe(+3) eaf788b fixed add_constant, undefined surface null pointer, added test cases 2212f9c fixed bug in reprep when sit had surface species. Added capability of sit + edl, have not tested it 79956e3 made tally_table a vector of class tally 58b0d1f Merge commit 'd77c11ec700085f19b76af6543013e23ee0739d3' d77c11e [phreeqci] fixed header error with phast 63175ab [phreeqci] fixed header error with phast 0feb715 [phreeqci] fixed WINDOWS.H already included error on windows builds 123cc8a [phreeqci] fixed _ASSERTE error on linux builds 22c4a62 [phreeqci] struct to class changes 4cee19d Merge commit '2d8ca2d0f37d13ad67be582208a4e65edfcf702f' 2d8ca2d [phreeqci] added 'new' debugging d0c8212 [phreeqci] added 'new' debugging 9661fea tokadd_heading leak 4565c5d catching upMerge remote-tracking branch 'origin/master' into classify c22d792 fix notab leak 6d2b45a Merge remote-tracking branch 'usgs-coupled/master' 38cfe18 memory leak user_print, pitz/sit store, add uphill_NPa, remove TESTINT 24f9bf7 removed TESTING definition e2ce928 Tony agreed with change for all_cells, new test case d2a5d63 reset all_cells in all cases e3c0d61 initialize aphi c960e05 builds on vs2005; still needs to initialize class pitz_param* aphi 71dc944 cl1mp, bad initialization 2e5f255 fixed errors/warnings from ming and intel 369733e converted to classes 7961b16 release.txt, couple size_t 5d76f82 copy operator works well enough 7ce8947 updated InternalCopy for operator equal 7bd13ff new/delete theta params, pitz_param_copy 50e8903 new/delete pitz_params 87d6792 reverting changes to sit_params and theta_params. Will consider using new and delet dcb9efe sit_params ac3335e theta_params 8878232 delete rate, unused cptr 492df61 descriptions 25e0621 cell_data 051ddba stag_data 33157a2 fixed more size_t and initialized all structs f86f430 back to original set of files I think af1b761 removing CReaction and Classes files 006d1de reorganizing 287f81c elt_list vectorized 7228bd0 move struct rxn_token 28de8b5 more size_t d2e3a4e Removed cxxChemRxn ce64720 cleaned up, removed struct reaction 028e908 moving to CReaction dc2dc53 vectorized token 9fd3f2a save_values rewritten with map 8a6cef5 vectorized save_values 8685225 fixed clang errors, needed .c_str 318e267 (size_t) max and count 1547d91 finished up spread b5c7ba4 going to work on warnings 4c848b4 all inverse structures vectorized. Starting on solver workspace 980d58e finished vectorizing struct inverse. Need to do sub structs d13bb76 removed count_elts 89ab28d vector inverse elts d575ade tidy.cpp, title_x 16fd18f removed string_duplicate from prep.cpp 82a10d6 revised get_elt and get_token d7e3be4 cleaned up some string_duplicate 76366a6 fixed processing file names 157a458 description_x 51fec19 class_main c748922 added const qualifier for all the parsing 380a6ea methods set to const, variables need to follow 6d67e22 copier and dash 48e6b93 fixed a new master, advection punch_temp and print_temp, some tidying 5f21daf unknown->master now a vector. Using size instead of a null to end list. 3c432d0 user_graph commands, alk_list 2b14f80 last_model 7a6b8b6 Merge branch 'warnings_redux' into vectorize_2 885a2f7 Fix memory bug in ex13_impl, tweak Makefile. 6907bb0 base, sit arrays 90e8412 starting on pitzer bd0cad9 vector kinetics arrays 1850c32 basic commands are now std::string 78a83ed c,d in polint d82d5d6 vector llnl parameters, removed hash references 7c538b6 Revert "delete s[i]" 97bcfd7 Merge branch 'warnings_redux' into vectorize 15a8991 delete s[i] 0b19404 master new/delete b100f85 more new/delete. Fixed str_tolower for ming fd93f84 needed to new/delete species and phase structs 1986e00 alphabetize tokens ee6fa53 bool analytic cc614e6 add_logk for logk, species, phases 67447c5 Removed hashtable, all hashes have been replaced with maps. ee7d2c5 replaced hash for isotope_ratio, isotope_alpha, calculate_value with maps. Fixed some case errors with new maps. 52e0622 replase master_isotope_hash_table with master_isotope_map c01c8d6 replace logk_hash_table with logk_map. Added str_tolower(std::string) 3e69461 replaced phases_hash_table with phases_map effafe0 replace species_hash_table with species_map 8bff6d3 removed HASH code. replaced elements_hash_table with elements_map 90e9ee0 removed ineq_init. Vector advection_print, advection_punch 2f38047 size_t for subscripts 5161ea7 Merged origin/master, Alphabetized Basic toks f8e05c1 only call qsort with more than one element 1ab8641 remove _v, use std::vector only, alloc at least 1 scratch 9732a1c cannot qsort size 0 vector 67fc478 one more .data 2f0f5e1 Some replacements of .data() were incorrect ba9813a remove .data() 43765f8 need <struct xxx> 0feb20d after merging origin/master, one fix needed f136feb Merging origin/master. Merge remote-tracking branch 'origin/master' into warnings_redux 71aa5b9 bug count_sys not incremented e43550c vector inverse d4cc14e vector x 6c0edef vector rates e3cc46a vector save_values 41b9965 vector species_list 449a54f vector mb_unknowns 51514eb vector delta, sum_jacobx f0707aa vector sum_mb1, sum_mb2 7d303de vector trxn.token 83cfb29 elt_list, moved qsort to elt_list_combine e8c9027 vector elt_list 0957a52 vector theta_params b1af156 vector pitz_params e3ea010 vector sit_params b87d0cd vector my_array, residual, delta e43471a vector s_x 622d361 vector s_x 3d41ef8 vector logk e8dd208 vector sys 3c9f594 vector master de1ba62 vector s e7c78a8 vector phases f2c64fe vector elements e8af689 vector isotope_alpha ba2601a vector isotope_ratio 76da4f8 finished master_isotope 4bb1c80 vector master_isotope 97e574d vector calculate_value** 9d9fbfb cl1 variables converted to std::vector 1e0d410 using memset 54b0d4d starting on space 5a649c2 Merge pull request #2 from usgs-coupled/gasphasepressures a992537 (void)sscanf, removed SKIP, removed PHREEQ98 6a5bb8a Merge pull request #1 from usgs-coupled/mar10 d9ced82 Fixed uninitialized constructors and couple of other warnings c79d2c2 working on UTF-8 fcee4d5 Added delta_h_species, delta_h_phase, dh_a0, dh_bdot Basic functions 81e862d Tonys changes Mar 10. SIs in inverse calulations 9e8b382 Merge remote-tracking branch 'usgs-coupled/master' 053b4c6 Merge remote-tracking branch 'origin/master' 20091aa Merge branch 'log10molalities' into gasphasepressures 41e1112 Last of changes for GetGasPhasePressures and GetGasPhasePhi, openmp and mpi. MPI fortrans not tested. e1f9cb1 more checking in. Should be down to tweaks for SetGasPhaseMoles. 00ee6e3 C++ is working with OpenMP and MPI for Get/SetGasPhaseMoles. Need to add c and F90. c3a3153 Added GetSpeciesLog10Molalities. Tested OpenMP with VS. Tested MPI with MinGW. Fortran, C, and C++ seem to work. e8b11f3 added optional 6th argument to Basic function sys to change sort order from molality or moles to the name. Added synonym PAD$. Added new mytest/sys_sort. 3e4fc7e cleanup commented lines 54b992f working on tabs and no newline 2181847 Merge branch 'master' of https://github.com/usgs-coupled/phreeqc3 deeecb0 needed strexpr in ADD_HEADING to allow expressions 9b7785f [iphreeqccom] updated date 711b1d0 Merge commit '608e74f5d3c55a4d91a4e08d86f2fd6df0ce0a05' 608e74f [wphast] updated date 5128e13 [phreeqc3] updated image location fba8ae2 [phreeqc3] updated image location 43988f0 initialize punch_newline 176fb02 Moved initialization from header to constructor, special characters in As.out c9f796a added ADD_HEADING for IPhreeqc 1362f0f Added EOL_NOTAB$ and NO_NEWLINE$, updated release notes 2b4dbbd Merge commit 'cd51d8aeed46909e5f028a19089acfef43d6ede9' f2023c4 Merge branch 'gtest' into 'master' cd51d8a reset for dlls 54161f4 reset for dlls 01c99a7 Merge remote-tracking branch 'github/master' 23f3917 Merge remote-tracking branch 'scharlton2/master' f6644e6 check for null pointer. Encoding for .out file 9319c9d Merge commit '5b816fa1fd82eb94e2702b6bd9df6066fb71267b' 5b816fa added src/phreeqcpp/common/PHRQ_exports.h 07717b1 added src/phreeqcpp/common/PHRQ_exports.h d8c638f Merge remote-tracking branch 'origin/master' into gtest 87bbb6a adjusted alignment for utf-8 strings 03bda16 added write_banner to non-DOS and added UTF8 define 995de52 converted to utf-8 fc8fe3e re-added src/ZedGraph.dll fbae3e9 code change for extending porosity definition. Change to TonyLitharge2a 46257e7 added googletest and fixed some minor bugs 13ca055 added googletest and fixed some minor bugs f1dda6c Fixed problem with exchange-related when exchanger is defined as CaX2 20daad4 I guess cxxSurface::NO_EDL is correct 801812d Tony's changes to implicit Nernst-Planck calculation 6b4892c added Basic function DEBYE_LENGTH and test case zeta 921ab10 Changed tidy_exchange_min and tidy_exchange_kin to tidy only for new_def and n_user >= 0. Fixed bug if surf_charge not defined for NO_EDL. Added test MoreExchMix 2aef60a Finished up surface and exchange related for cases where related phase or kinetic reaction was modified. Proportionality should now be maintained. Added test cases. 569e1e1 Exchange related. Needed to update in case the related entity changed. ea54e02 Free str in callback in PBasic a87cd1f Merge commit '1871b026ca8487c23a025415dbc0b2eca01f9af4' 1871b02 fixed some c2011 warnings, added more info for -formula errors, fixed pressure llnl examples aa4d023 fixed some c2011 warnings, added more info for -formula errors, fixed pressure llnl examples e1465e3 Commit from David's Email 2020-05-22; Implementing llnl-type databases with higher temperature nad pressure e18e1ec Tony bug fix for TRANSPORT. Harmonic mean for boundary? Added Cub example. 44f077e Merge commit 'e68934133fc9cd45e7cccc397c55e13f7ee92e5b' e689341 [phreeqci] Testing subtree merges 4f34fd0 [phreeqci] Testing subtree merges 69c0bb3 fixed conflict on merge 55c4dba Merge commit 'b25fc5bdd48b6d3ab8d677f7d38ad3a462789500' b25fc5b fixed conflict on merge ca80be6 fixed conflict on merge 49a74a6 [phreeqc3] Testing subtree merges aec6f90 [phreeqc3] Testing subtree merges c4c224a Merge commit '84865ad5ac30a9edb86c89ced4194d127ee896fd' 0bf4138 Merge commit '4a8727cecd9fefd1587485820e913c0e666b77d9' 553875f Merge commit 'aab8bc12ea8be8aec5943e1c77a54b19d28168cb' aab8bc1 Merge commit '84865ad5ac30a9edb86c89ced4194d127ee896fd' 7bd02ff Fixed bug with more porosities than cells in TRANSPORT. Added silica sorption to databases. Revised CalPortDiff 84865ad Added .gitlab-ci.yml d398195 Added .gitlab-ci.yml 40c2787 Added .gitlab-ci.yml 3b6ce6c Added .gitlab-ci.yml daf64a1 Added .gitlab-ci.yml ae06f35 Fixed GFW bug on new elements in TRANSPORT 9cc783b added Basic functions for PHAST: velocity_x, velocity_y, velocity_z, transport_cell_no 79f768a Merge branch 'master' into 'master' bd7634a removed j = j in loop 542394c IPhreeqc: ifdef'd out references to std::cerr and std::cout 6067ce8 Merge branch 'implicit3' into 'master' 21bd20f Fixed more compile warnings. Removed andra_kin_ss from testing, results are inconsistent between Linux and Windows, presumably the ifs in RATES 97b9c58 Merge branch 'implicit3' into 'master' 45db5cf Another Linux warning, lower tol on andra_kin_ss. 443be1c Merge branch 'implicit3' into 'master' 9a29aaf Last Linux compile warnings. Added more precision to andra_kin_ss. 6dafd7d Merge branch 'implicit3' into 'master' fbde633 Fixing Linux compiler warnings, checking in new regression test files. 2207711 Merge branch 'implicit3' into 'master' 77e36a2 Tony fixed some transport, revised colloid_U. New cases added to Makefile. f07caf9 Changed back print to allow incremental_reactions to work correctly beadd07 Merge commit '5947da90657d1cb8f832152b4573dca0bbefb49e' 6a49d41 changes to make related and mixing items independent of case. surface_mix test case. 5947da9 initial Tony changes 8089c10 initial Tony changes 009aec7 Merge remote-tracking branch 'coupled/master' 4676ee4 added more P-R gas paramteters c07314c Merge commit '492a4d257f300b7a9e0b5dc7e212c8f85ecb7f6e' 492a4d2 Merge remote-tracking branch 'coupled/master' 81ca633 Merge remote-tracking branch 'coupled/master' 950fca2 CRAN: replaced deprecated std::ptr_fun with lambda function 597bcd7 CRAN: replaced deprecated std::ptr_fun with lambda function 044e0ea phreeqc_ptr bug in internal copy 5934297 Merge commit '5c53fb207238bc0e846123a7e0d71a48bd9976ab' 5c53fb2 Merge commit '1327e93127e40e7a55ec629dcc9dd91ec29e77fe' c117e18 Tony fix of index error b90ddb5 Fixed Tony's fix, added implicit_as example 03acc3f changed abs to fabs 1fef40e added implicit, max_mixf to internal copy 32939ba Merge commit '1327e93127e40e7a55ec629dcc9dd91ec29e77fe' b3bf691 fixed > > in templates for gcc c929113 Tony fix May 31 1327e93 Implicit seems to be working with Tony's latest changes 55ea163 Implicit seems to be working with Tony's latest changes c7111f7 Sort of works, still bugs and serious errors compared to explicit 600c7ee Fixed some bugs with iso.dat inverse modeling, added test case. Still does not generate [13C](4) and [13C](-4) from SOLUTION 2291700 Fixed gas_phase_mix bug, added test case 035a4e0 Tony tweak to transport.cpp bd4fc25 Merge branch 'tony20190117' into 'master' 71c994b skipping restart 1257f8c Merge branch 'issue-3' into 'master' ce33478 Fixed -Wcatch-value warnings reported by CRAN 040fd95 include restart, remove ex20_debug d57264d 2. changes to solid solution numerical method 3fd8155 changes to solid solution numerical method 2b14a94 Tony's changes 20190117 ae6e8b0 added modify methods for restart files b500c54 changed restart file to include UZ fffac6d another try for ex20_debug fa5ee50 fix problem with ex20_debug d993901 encoding, limit.h 92c81f9 Revised logic for nmix 3cc84da Merge remote-tracking branch 'coupled/master' into ss_trans 56b5bf3 create valid ranges when sscanf doesn't return 2 c43c9af tweaked ss, changed surf function per Kinniburgh b10df16 Corrected syntax of integer limit, previous commit actually changed ss convergence parameter, used to multiply by 0.99 d74c8ff Corrected syntax of integer limit 906cfd4 Check value of nmix 058375c removed check of ss when sum of components is small 2977db4 Tonys fix to diffusion bug with porosity change f904467 revised lists to be cumulative for eq, gas, kin, ss 9285985 merging coupled/master into copy 7c23b62 Fixed string_duplicate memory error 2d5551a fixed sc7 for copy and initial time 4842d9e inverse iter 100000; finished copy operator; a bit more testing to go 4eefe43 ex20_debug fix 78e39cd still debugging copy cee10e7 fixing bugs in copy operator ebab4bc fixing bugs in copy operator 5a35e02 Fixed Linux warnings, memory errors b86f793 Beginning to test copy operator 5d40e69 [IPhreeqc] added parens for clang++ -Wlogical-op-parentheses 936de38 removed register keywords and updated for misc clang warnings ec9de4c beginning of checking copy operation ebeddcd [iphreeqc] Changes for CRAN 3.4.7 9592d6e Merge branch 'dlpark-phreeqc3-TonyApr2018' into 'master' 7c0fb65 [phreeqc3] needed to check gas phase type for same model, added test case 9152ca2 Closes #1 ebc4f69 Merge branch 'dlpark-phreeqc3-TonyApr2018' into 'master' 97a0cec Fixed bug where 1W was interpreted as an isotope 2deb4ed added option -ddl to surface. Added test case df7d5de Merge branch 'gammas' into 'master' 34abb5b gammas finished, working on reactants 5314827 Tony's changes; diffuse layer with pitzer 4271ca4 Tonys corrections, added balonis test 2e390fd commit fix for Mtg git-subtree-dir: phreeqcpp git-subtree-split: 2243d25babbc524e7875b3d591bb6b91c4399a95
5154 lines
136 KiB
C++
5154 lines
136 KiB
C++
#include "Phreeqc.h"
|
|
#include "phqalloc.h"
|
|
#include "Utils.h"
|
|
#include "Solution.h"
|
|
#include "SolutionIsotope.h"
|
|
|
|
|
|
#define MAX_MODELS 20
|
|
#define MIN_TOTAL_INVERSE 1e-14
|
|
|
|
/* variables local to module */
|
|
#define SCALE_EPSILON .0009765625
|
|
#define SCALE_WATER 1.
|
|
#define SCALE_ALL 1.
|
|
|
|
#if defined(PHREEQCI_GUI)
|
|
#ifdef _DEBUG
|
|
#define new DEBUG_NEW
|
|
#undef THIS_FILE
|
|
static char THIS_FILE[] = __FILE__;
|
|
#endif
|
|
#endif
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
inverse_models(void)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Go through list of inverse models, make calculations
|
|
* for any marked "new".
|
|
*/
|
|
int n/*, print1*/;
|
|
char string[MAX_LENGTH];
|
|
if (count_inverse <= 0) return OK;
|
|
// Revert to previous headings after inverse modeling
|
|
std::vector<std::string> old_headings;
|
|
|
|
state = INVERSE;
|
|
dl_type_x = cxxSurface::NO_DL;
|
|
|
|
for (n = 0; n < count_inverse; n++)
|
|
{
|
|
if (inverse[n].new_def == TRUE)
|
|
{
|
|
/*
|
|
* dump .lon file
|
|
*/
|
|
if (inverse[n].netpath != NULL)
|
|
dump_netpath(&inverse[n]);
|
|
|
|
/*
|
|
* open .pat file
|
|
*/
|
|
if (inverse[n].pat != NULL)
|
|
{
|
|
strcpy(string, inverse[n].pat);
|
|
if (replace(".pat", ".pat", string) != TRUE)
|
|
{
|
|
strcat(string, ".pat");
|
|
}
|
|
netpath_file = fopen(string, "w");
|
|
if (netpath_file == NULL)
|
|
{
|
|
error_string = sformatf( "Can`t open file, %s.", string);
|
|
error_msg(error_string, STOP);
|
|
#if !defined(R_SO)
|
|
exit(4);
|
|
#endif
|
|
}
|
|
count_inverse_models = 0;
|
|
count_pat_solutions = 0;
|
|
/* Header */
|
|
fprintf(netpath_file, "2.14 # File format\n");
|
|
}
|
|
/*
|
|
* Fill in stucture "use".
|
|
*/
|
|
use.Set_inverse_in(true);
|
|
use.Set_inverse_ptr(&inverse[n]);
|
|
use.Set_n_inverse_user(inverse[n].n_user);
|
|
/*
|
|
* Initial prints
|
|
*/
|
|
error_string = sformatf(
|
|
"Beginning of inverse modeling %d calculations.",
|
|
inverse[n].n_user);
|
|
dup_print(error_string, TRUE);
|
|
|
|
if (inverse[n].mp == TRUE)
|
|
{
|
|
output_msg(sformatf(
|
|
"Using Cl1MP multiprecision optimization routine.\n"));
|
|
|
|
}
|
|
else
|
|
{
|
|
output_msg(sformatf(
|
|
"Using Cl1 standard precision optimization routine.\n"));
|
|
}
|
|
status(0, NULL);
|
|
|
|
/*
|
|
* Setup and solve
|
|
*/
|
|
count_calls = 0;
|
|
setup_inverse(&(inverse[n]));
|
|
punch_model_heading(&inverse[n]);
|
|
solve_inverse(&(inverse[n]));
|
|
if (inverse[n].isotope_unknowns.size() > 0)
|
|
{
|
|
inverse[n].isotope_unknowns.clear();
|
|
}
|
|
inverse[n].new_def = FALSE;
|
|
if (inverse[n].pat != NULL)
|
|
{
|
|
fclose(netpath_file);
|
|
netpath_file = NULL;
|
|
}
|
|
}
|
|
}
|
|
|
|
//user_punch_count_headings = (int) old_headings.size();
|
|
//user_punch_headings = (const char **) PHRQ_realloc(user_punch_headings,
|
|
// (size_t) (user_punch_count_headings + 1) * sizeof(char *));
|
|
//if (user_punch_headings == NULL)
|
|
// malloc_error();
|
|
//for (i = 0; i < user_punch_count_headings; i++)
|
|
//{
|
|
// user_punch_headings[i] = string_hsave(old_headings[i].c_str());
|
|
//}
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
setup_inverse(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Fill in array for an inverse problem
|
|
*/
|
|
int i, j, k, i_alk, i_carb;
|
|
size_t max;
|
|
size_t count_rows_t;
|
|
size_t column, row;
|
|
int temp;
|
|
LDBLE isotope_number;
|
|
LDBLE f, coef, cb, conc;
|
|
char token[MAX_LENGTH];
|
|
class phase *phase_ptr;
|
|
cxxSolution *solution_ptr;
|
|
CReaction *rxn_ptr;
|
|
class master *master_ptr;
|
|
/*
|
|
* Determine array sizes, row and column positions
|
|
*/
|
|
toler = inv_ptr->tolerance;
|
|
if (inv_ptr->mp == TRUE)
|
|
{
|
|
toler = inv_ptr->mp_tolerance;
|
|
}
|
|
/*
|
|
* Alkalinity derivatives with pH and carbon
|
|
*/
|
|
carbon = 1;
|
|
temp = pr.status;
|
|
pr.status = FALSE;
|
|
|
|
// current_selected_output is NULL at this point
|
|
carbon_derivs(inv_ptr);
|
|
pr.status = temp;
|
|
|
|
//current_selected_output->Set_inverse(temp_inv);
|
|
state = INVERSE;
|
|
/*
|
|
* tidy isotopes if necessary
|
|
*/
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
set_isotope_unknowns(inv_ptr);
|
|
if (get_input_errors() > 0)
|
|
{
|
|
error_msg("Stopping because of input errors.", STOP);
|
|
}
|
|
check_isotopes(inv_ptr);
|
|
if (get_input_errors() > 0)
|
|
{
|
|
error_msg("Stopping because of input errors.", STOP);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* count unknowns
|
|
*/
|
|
max_column_count = inv_ptr->elts.size() * inv_ptr->count_solns + /* epsilons */
|
|
inv_ptr->count_solns + /* solutions */
|
|
inv_ptr->phases.size() + /* phases */
|
|
inv_ptr->count_redox_rxns + /* redox reactions */
|
|
carbon * inv_ptr->count_solns + /* pH */
|
|
1 + /* water */
|
|
inv_ptr->isotope_unknowns.size() * inv_ptr->count_solns + /* isotopes in solution */
|
|
inv_ptr->isotopes.size() * inv_ptr->phases.size() + /* isotopes in phases */
|
|
1 + 1; /* rhs, ineq */
|
|
count_unknowns = max_column_count - 2;
|
|
col_phases = inv_ptr->count_solns;
|
|
col_redox = col_phases + inv_ptr->phases.size();
|
|
col_epsilon = col_redox + inv_ptr->count_redox_rxns;
|
|
col_ph = col_epsilon + inv_ptr->elts.size() * inv_ptr->count_solns;
|
|
col_water = col_ph + carbon * inv_ptr->count_solns;
|
|
col_isotopes = col_water + 1;
|
|
col_phase_isotopes =
|
|
col_isotopes + inv_ptr->isotope_unknowns.size() * inv_ptr->count_solns;
|
|
max_row_count = inv_ptr->count_solns * inv_ptr->elts.size() + /* optimize */
|
|
carbon * inv_ptr->count_solns + /* optimize ph */
|
|
1 + /* optimize water */
|
|
inv_ptr->count_solns * inv_ptr->isotope_unknowns.size() + /* optimize isotopes */
|
|
inv_ptr->isotopes.size() * inv_ptr->phases.size() + /* optimize phase isotopes */
|
|
inv_ptr->elts.size() + /* mass balances */
|
|
1 + 1 + /* fractions, init and final */
|
|
inv_ptr->count_solns + /* charge balances */
|
|
carbon * inv_ptr->count_solns + /* dAlk = dC + dph */
|
|
inv_ptr->isotopes.size() + /* isotopes */
|
|
2 * inv_ptr->count_solns * inv_ptr->elts.size() + /* epsilon constraints */
|
|
2 * carbon * inv_ptr->count_solns + /* epsilon on ph */
|
|
2 + /* epsilon for water */
|
|
2 * inv_ptr->isotope_unknowns.size() * inv_ptr->count_solns + /* epsilon for isotopes */
|
|
2 * inv_ptr->isotopes.size() * inv_ptr->phases.size() + /* epsilon for isotopes in phases */
|
|
2; /* work space */
|
|
|
|
row_mb = inv_ptr->count_solns * inv_ptr->elts.size() +
|
|
carbon * inv_ptr->count_solns + 1 +
|
|
inv_ptr->count_solns * inv_ptr->isotope_unknowns.size() +
|
|
inv_ptr->isotopes.size() * inv_ptr->phases.size();
|
|
row_fract = row_mb + inv_ptr->elts.size();
|
|
row_charge = row_fract + 2;
|
|
row_carbon = row_charge + inv_ptr->count_solns;
|
|
row_isotopes = row_carbon + carbon * inv_ptr->count_solns;
|
|
row_epsilon = row_isotopes + inv_ptr->isotopes.size();
|
|
/* The next three are not right, some rows of epsilon are deleted */
|
|
/*
|
|
row_ph_epsilon = row_epsilon + 2 * inv_ptr->count_solns * inv_ptr->count_elts;
|
|
row_water_epsilon = row_ph + 2 * carbon * inv_ptr->count_solns;
|
|
row_isotope_epsilon
|
|
row_isotope_phase_epsilon
|
|
*/
|
|
/*
|
|
* Malloc space for arrays
|
|
*/
|
|
my_array.resize(max_column_count * max_row_count);
|
|
array1.resize(max_column_count * max_row_count);
|
|
col_name.resize(max_column_count);
|
|
row_name.resize(max_row_count);
|
|
delta.resize(max_column_count);
|
|
inv_delta1.resize(max_column_count);
|
|
delta2.resize(max_column_count);
|
|
delta3.resize(max_column_count);
|
|
delta_save.resize(max_column_count);
|
|
min_delta.resize(max_column_count);
|
|
max_delta.resize(max_column_count);
|
|
inv_res.resize(max_row_count);
|
|
|
|
if (max_column_count < max_row_count)
|
|
{
|
|
max = max_row_count;
|
|
}
|
|
else
|
|
{
|
|
max = max_column_count;
|
|
}
|
|
inv_zero.resize((size_t) max);
|
|
/*
|
|
* Define inv_zero and inv_zero array, delta
|
|
*/
|
|
for (i = 0; i < (int) max; i++)
|
|
inv_zero[i] = 0.0;
|
|
|
|
memcpy((void *) &(delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(min_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(max_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
for (i = 0; i < max_row_count; i++)
|
|
{
|
|
memcpy((void *) &(my_array[(size_t)i * max_column_count]), (void *) &(inv_zero[0]),
|
|
max_column_count * sizeof(LDBLE));
|
|
}
|
|
/*
|
|
* begin filling array
|
|
*/
|
|
count_rows = 0;
|
|
/*
|
|
* optimization
|
|
*/
|
|
count_optimize = inv_ptr->count_solns * inv_ptr->elts.size() + /* optimize */
|
|
carbon * inv_ptr->count_solns + /* optimize ph */
|
|
1 + /* optimize water */
|
|
inv_ptr->count_solns * inv_ptr->isotope_unknowns.size() + /* optimize isotopes */
|
|
inv_ptr->isotopes.size() * inv_ptr->phases.size(); /* optimize phase isotopes */
|
|
|
|
for (i = 0; i < count_optimize; i++)
|
|
{
|
|
row_name[count_rows] = string_hsave("optimize");
|
|
count_rows++;
|
|
}
|
|
write_optimize_names(inv_ptr);
|
|
/*
|
|
* equalities
|
|
*/
|
|
|
|
/*
|
|
* Mass_balance: solution data
|
|
*/
|
|
|
|
/* initialize master species */
|
|
for (i = 0; i < (int)master.size(); i++)
|
|
{
|
|
master[i]->in = -1;
|
|
if (strstr(master[i]->elt->name, "Alk") == master[i]->elt->name)
|
|
{
|
|
master_alk = master[i];
|
|
}
|
|
}
|
|
/* mark master species included in model, write row names */
|
|
count_rows_t = count_rows;
|
|
i_alk = -1;
|
|
i_carb = -1;
|
|
for (i = 0; i < inv_ptr->elts.size(); i++)
|
|
{
|
|
master_ptr = inv_ptr->elts[i].master;
|
|
if (master_ptr == master_alk)
|
|
i_alk = i;
|
|
if (strcmp(master_ptr->elt->name, "C(4)") == 0)
|
|
i_carb = i;
|
|
inv_ptr->elts[i].master->in = (int)count_rows_t;
|
|
row_name[count_rows_t] = inv_ptr->elts[i].master->elt->name;
|
|
count_rows_t++;
|
|
}
|
|
/* put concentrations in array */
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
xsolution_zero();
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
if (solution_ptr == NULL)
|
|
{
|
|
error_string = sformatf( "Solution number %d not found.",
|
|
inv_ptr->solns[i]);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
/* write master species concentrations */
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
master_ptr = master_bsearch(jit->first.c_str());
|
|
master_ptr->total += jit->second;
|
|
/* List elements not included in model */
|
|
if (master_ptr->in < 0)
|
|
{
|
|
error_string = sformatf(
|
|
"%s is included in solution %d, but is not included as a mass-balance constraint.",
|
|
jit->first.c_str(),
|
|
inv_ptr->solns[i]);
|
|
warning_msg(error_string);
|
|
}
|
|
}
|
|
master_alk->total = solution_ptr->Get_total_alkalinity();
|
|
f = 1.0;
|
|
if (i == (inv_ptr->count_solns - 1))
|
|
{
|
|
f = -1.0;
|
|
}
|
|
column = i;
|
|
sprintf(token, "soln %d", i);
|
|
col_name[column] = string_hsave(token);
|
|
for (j = 0; j < (int)master.size(); j++)
|
|
{
|
|
if (master[j]->in >= 0)
|
|
{
|
|
my_array[(size_t)master[j]->in * max_column_count + (size_t)i] =
|
|
f * master[j]->total;
|
|
if (master[j]->s == s_eminus)
|
|
{
|
|
my_array[(size_t)master[j]->in * max_column_count + (size_t)i] = 0.0;
|
|
}
|
|
}
|
|
}
|
|
/* calculate charge balance for elements in model */
|
|
cb = 0;
|
|
for (j = 0; j < (int)master.size(); j++)
|
|
{
|
|
if (master[j]->in >= 0)
|
|
{
|
|
if (master[j]->s == s_eminus)
|
|
{
|
|
coef = 0.0;
|
|
}
|
|
else if (master[j] == master_alk)
|
|
{
|
|
coef = -1.0;
|
|
}
|
|
else
|
|
{
|
|
coef = master[j]->s->z + master[j]->s->alk;
|
|
}
|
|
cb += coef * master[j]->total;
|
|
}
|
|
}
|
|
if (fabs(cb) < toler)
|
|
cb = 0.0;
|
|
my_array[((size_t)row_charge + (size_t)i) * max_column_count + (size_t)i] = cb;
|
|
}
|
|
|
|
/* mass_balance: phase data */
|
|
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
phase_ptr = inv_ptr->phases[i].phase;
|
|
rxn_ptr = &phase_ptr->rxn_s;
|
|
column = col_phases + i;
|
|
col_name[column] = phase_ptr->name;
|
|
for (j = 1; rxn_ptr->token[j].s != NULL; j++)
|
|
{
|
|
if (rxn_ptr->token[j].s->secondary != NULL)
|
|
{
|
|
master_ptr = rxn_ptr->token[j].s->secondary;
|
|
}
|
|
else
|
|
{
|
|
master_ptr = rxn_ptr->token[j].s->primary;
|
|
}
|
|
if (master_ptr == NULL)
|
|
{
|
|
error_string = sformatf(
|
|
"Setup_inverse, reaction for phase, %s.",
|
|
phase_ptr->name);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
if (master_ptr->s == s_hplus)
|
|
continue;
|
|
if (master_ptr->s == s_h2o)
|
|
{
|
|
row = row_fract;
|
|
/* turn off h2o from minerals in water mass balance */
|
|
if (inv_ptr->mineral_water != TRUE)
|
|
continue;
|
|
|
|
}
|
|
else
|
|
{
|
|
row = master_ptr->in;
|
|
}
|
|
/* e- has coef of 0 for some reason */
|
|
coef = master_ptr->coef;
|
|
if (coef <= 0)
|
|
coef = 1.0;
|
|
my_array[(size_t)row * max_column_count + (size_t)column] =
|
|
rxn_ptr->token[j].coef * coef;
|
|
}
|
|
row = master_alk->in; /* include alkalinity for phase */
|
|
my_array[(size_t)row * max_column_count + (size_t)column] = calc_alk(*rxn_ptr);
|
|
}
|
|
|
|
/* mass balance: redox reaction data */
|
|
|
|
k = 0;
|
|
for (i = 0; i < inv_ptr->elts.size(); i++)
|
|
{
|
|
if (inv_ptr->elts[i].master->s->primary == NULL)
|
|
{
|
|
coef = inv_ptr->elts[i].master->coef;
|
|
rxn_ptr = &inv_ptr->elts[i].master->rxn_primary;
|
|
column = col_redox + k;
|
|
col_name[column] = inv_ptr->elts[i].master->elt->name;
|
|
k++;
|
|
for (j = 0; rxn_ptr->token[j].s != NULL; j++)
|
|
{
|
|
if (rxn_ptr->token[j].s->secondary != NULL)
|
|
{
|
|
master_ptr = rxn_ptr->token[j].s->secondary;
|
|
}
|
|
else
|
|
{
|
|
master_ptr = rxn_ptr->token[j].s->primary;
|
|
}
|
|
if (master_ptr == NULL)
|
|
{
|
|
error_string = sformatf(
|
|
"Subroutine setup_inverse, element not found, %s.",
|
|
rxn_ptr->token[j].s->name);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
if (master_ptr->s == s_hplus)
|
|
continue;
|
|
if (master_ptr->s == s_h2o)
|
|
{
|
|
row = row_fract;
|
|
/* turn off h2o from minerals in water mass balance */
|
|
if (inv_ptr->mineral_water != TRUE)
|
|
continue;
|
|
}
|
|
else
|
|
{
|
|
row = master_ptr->in;
|
|
}
|
|
assert(row * max_column_count + column < max_column_count * max_row_count);
|
|
assert(row >= 0);
|
|
assert(column >= 0);
|
|
my_array[(size_t)row * max_column_count + (size_t)column] =
|
|
rxn_ptr->token[j].coef;
|
|
/* if coefficient of element is not 1.0 in master species */
|
|
if (j != 0)
|
|
my_array[(size_t)row * max_column_count + (size_t)column] /= coef;
|
|
}
|
|
row = master_alk->in; /* include alkalinity for redox reaction */
|
|
my_array[(size_t)row * max_column_count + (size_t)column] =
|
|
(calc_alk(*rxn_ptr) - inv_ptr->elts[i].master->s->alk) / coef;
|
|
}
|
|
}
|
|
|
|
/* mass-balance: epsilons */
|
|
|
|
column = col_epsilon;
|
|
for (i = 0; i < inv_ptr->elts.size(); i++)
|
|
{
|
|
row = inv_ptr->elts[i].master->in;
|
|
for (j = 0; j < inv_ptr->count_solns; j++)
|
|
{
|
|
if (j < (inv_ptr->count_solns - 1))
|
|
{
|
|
my_array[(size_t)row * max_column_count + (size_t)column] = 1.0;
|
|
}
|
|
else
|
|
{
|
|
my_array[(size_t)row * max_column_count + (size_t)column] = -1.0;
|
|
}
|
|
if (inv_ptr->elts[i].master->s == s_eminus)
|
|
{
|
|
my_array[(size_t)row * max_column_count + (size_t)column] = 0.0;
|
|
}
|
|
sprintf(token, "%s %d", row_name[row], j);
|
|
col_name[column] = string_hsave(token);
|
|
column++;
|
|
}
|
|
}
|
|
count_rows += inv_ptr->elts.size();
|
|
|
|
/* put names in col_name for ph */
|
|
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
sprintf(token, "ph %d", i);
|
|
col_name[column] = string_hsave(token);
|
|
column++;
|
|
}
|
|
/* put names in col_name for water */
|
|
|
|
sprintf(token, "water");
|
|
col_name[column] = string_hsave(token);
|
|
column++;
|
|
|
|
/* put names of isotopes in col_name */
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
for (j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
sprintf(token, "%d%s %d",
|
|
(int) inv_ptr->isotope_unknowns[j].isotope_number,
|
|
inv_ptr->isotope_unknowns[j].elt_name, i);
|
|
col_name[column] = string_hsave(token);
|
|
column++;
|
|
}
|
|
}
|
|
|
|
/* put phase isotopes in col_name */
|
|
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
/* isotopes of phases phases */
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
for (j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
sprintf(token, "%d%s %s",
|
|
(int) inv_ptr->isotopes[j].isotope_number,
|
|
inv_ptr->isotopes[j].elt_name,
|
|
inv_ptr->phases[i].phase->name);
|
|
col_name[column] = string_hsave(token);
|
|
column++;
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Initial solution mixing fractions or water mass balance
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
if (i < inv_ptr->count_solns - 1)
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)i] =
|
|
1.0 / gfw_water * solution_ptr->Get_mass_water();
|
|
}
|
|
else
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)inv_ptr->count_solns - 1] =
|
|
-1.0 / gfw_water * solution_ptr->Get_mass_water();
|
|
}
|
|
}
|
|
/* coefficient for water uncertainty */
|
|
if (inv_ptr->water_uncertainty > 0)
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)col_water] = 1.0;
|
|
}
|
|
row_name[count_rows] = string_hsave("H2O");
|
|
row_water = count_rows;
|
|
count_rows++;
|
|
|
|
/*
|
|
* Final solution fraction equals 1.0
|
|
*/
|
|
|
|
my_array[count_rows * max_column_count + (size_t)inv_ptr->count_solns - 1] = 1.0;
|
|
my_array[count_rows * max_column_count + count_unknowns] = 1.0;
|
|
row_name[count_rows] = string_hsave("fract, final");
|
|
count_rows++;
|
|
|
|
/*
|
|
* Charge balance:
|
|
*/
|
|
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
/* solution_ptr = solution_bsearch(inv_ptr->solns[i], &j, TRUE); */
|
|
/* array[count_rows * max_column_count + (size_t)i] = solution_ptr->cb; */
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
column = col_epsilon + j * inv_ptr->count_solns + i;
|
|
coef =
|
|
inv_ptr->elts[j].master->s->z +
|
|
inv_ptr->elts[j].master->s->alk;
|
|
if (inv_ptr->elts[j].master == master_alk)
|
|
{
|
|
coef = -1.0;
|
|
}
|
|
my_array[count_rows * max_column_count + (size_t)column] = coef;
|
|
if (inv_ptr->elts[j].master->s == s_eminus)
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)column] = 0.0;
|
|
}
|
|
}
|
|
sprintf(token, "%s %d", "charge", i);
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
/*
|
|
* dC = (dC/dph)*dph + (dC/dAlk)*dAlk for each solution
|
|
*/
|
|
/*
|
|
* dAlk = (dAlk/dC)*dC + (dAlk/dpH)*dpH for each solution
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (inv_ptr->dalk_dph[i] != 0 || inv_ptr->dalk_dc[i] != 0)
|
|
{
|
|
column = col_ph + i;
|
|
my_array[count_rows * max_column_count + (size_t)column] =
|
|
inv_ptr->dalk_dph[i];
|
|
column = col_epsilon + i_alk * inv_ptr->count_solns + i;
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
column = col_epsilon + i_carb * inv_ptr->count_solns + i;
|
|
my_array[count_rows * max_column_count + (size_t)column] =
|
|
inv_ptr->dalk_dc[i];
|
|
}
|
|
sprintf(token, "%s %d", "dAlk", i);
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
/*
|
|
* Isotope mass balances
|
|
*/
|
|
if (get_input_errors() > 0)
|
|
{
|
|
error_msg("Stopping because of input errors.", STOP);
|
|
}
|
|
if (inv_ptr->isotopes.size() != 0)
|
|
{
|
|
for (size_t j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
isotope_balance_equation(inv_ptr, (int)count_rows, (int)j);
|
|
sprintf(token, "%d%s", (int) inv_ptr->isotopes[j].isotope_number,
|
|
inv_ptr->isotopes[j].elt_name);
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
}
|
|
/*
|
|
* inequalities
|
|
*/
|
|
row_epsilon = count_rows;
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
if (inv_ptr->elts[j].master->s == s_eminus)
|
|
continue;
|
|
column = col_epsilon + j * inv_ptr->count_solns + i;
|
|
|
|
/* calculate magnitude of bound */
|
|
|
|
coef = inv_ptr->elts[j].uncertainties[i];
|
|
if (coef <= 0.0)
|
|
{
|
|
coef = -coef;
|
|
}
|
|
else
|
|
{
|
|
coef = my_array[(size_t)inv_ptr->elts[j].master->in *
|
|
max_column_count + (size_t)i] * coef;
|
|
coef = fabs(coef);
|
|
}
|
|
|
|
if (coef < toler)
|
|
coef = 0;
|
|
|
|
/* zero column if uncertainty is zero */
|
|
if (coef == 0.0)
|
|
{
|
|
for (k = 0; k < count_rows; k++)
|
|
{
|
|
my_array[(size_t)k * max_column_count + (size_t)column] = 0.0;
|
|
}
|
|
continue;
|
|
}
|
|
|
|
/* this statement probably obviates some of the following logic. */
|
|
/* coef += toler; */
|
|
|
|
/* scale epsilon optimization equation */
|
|
|
|
if (coef < toler)
|
|
{
|
|
my_array[((size_t)column - (size_t)col_epsilon) * max_column_count + (size_t)column] =
|
|
SCALE_EPSILON / toler;
|
|
}
|
|
else
|
|
{
|
|
my_array[((size_t)column - (size_t)col_epsilon) * max_column_count + (size_t)column] =
|
|
SCALE_EPSILON / coef;
|
|
}
|
|
|
|
/* set upper limit of change in positive direction */
|
|
if (coef < toler)
|
|
{
|
|
coef = toler;
|
|
f = 10;
|
|
}
|
|
else
|
|
{
|
|
f = 1.0;
|
|
}
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0 * f;
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef * f;
|
|
sprintf(token, "%s %s", inv_ptr->elts[j].master->elt->name, "eps+");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* set lower limit of change in negative direction */
|
|
conc = my_array[(size_t)inv_ptr->elts[j].master->in * max_column_count + (size_t)i];
|
|
|
|
/* if concentration is zero, only positive direction allowed */
|
|
if (conc == 0.0)
|
|
{
|
|
delta[column] = 1.0;
|
|
continue;
|
|
}
|
|
/* if uncertainty is less than tolerance, set uncertainty to toler */
|
|
if (coef <= toler)
|
|
{
|
|
/* f = 10 * toler / coef; */
|
|
coef = toler;
|
|
f = 10;
|
|
}
|
|
else
|
|
{
|
|
f = 1.0;
|
|
}
|
|
/* if uncertainty is greater than concentration,
|
|
maximum negative is equal to concentrations,
|
|
except alkalinity */
|
|
if (coef > fabs(conc) &&
|
|
(strstr(inv_ptr->elts[j].master->elt->name, "Alkalinity") !=
|
|
inv_ptr->elts[j].master->elt->name))
|
|
coef = fabs(conc) + toler;
|
|
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef * f;
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0 * f;
|
|
sprintf(token, "%s %s", inv_ptr->elts[j].master->elt->name,
|
|
"eps-");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
}
|
|
/*
|
|
* inequalities for pH
|
|
*/
|
|
/* row_ph_epsilon = count_rows; */
|
|
if (inv_ptr->carbon == TRUE)
|
|
{
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
column = col_ph + i;
|
|
coef = inv_ptr->ph_uncertainties[i];
|
|
|
|
/* scale epsilon in optimization equation */
|
|
|
|
my_array[((size_t)column - (size_t)col_epsilon) * max_column_count + (size_t)column] =
|
|
SCALE_EPSILON / coef;
|
|
|
|
/* set upper limit of change in positive direction */
|
|
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0;
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef;
|
|
sprintf(token, "%s %s", "pH", "eps+");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* set lower limit of change in negative direction */
|
|
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef;
|
|
sprintf(token, "%s %s", "pH", "eps-");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
}
|
|
/*
|
|
* inequalities for water
|
|
*/
|
|
column = col_water;
|
|
coef = inv_ptr->water_uncertainty;
|
|
/* row_water_epsilon = count_rows; */
|
|
if (coef > 0.0)
|
|
{
|
|
/* set upper limit of change in positive direction */
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0;
|
|
my_array[count_rows * max_column_count + count_unknowns] = coef;
|
|
sprintf(token, "%s %s", "water", "eps+");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* set lower limit of change in negative direction */
|
|
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
my_array[count_rows * max_column_count + count_unknowns] = coef;
|
|
sprintf(token, "%s %s", "water", "eps-");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
}
|
|
/*
|
|
* inequalities for isotopes
|
|
*/
|
|
row_isotope_epsilon = count_rows;
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
for (j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
column =
|
|
col_isotopes + (i * inv_ptr->isotope_unknowns.size()) + j;
|
|
master_ptr = inv_ptr->isotope_unknowns[j].master;
|
|
isotope_number = inv_ptr->isotope_unknowns[j].isotope_number;
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
class master *master_kit = master_bsearch(kit->second.Get_elt_name().c_str());
|
|
if (master_kit == master_ptr &&
|
|
kit->second.Get_isotope_number() ==
|
|
isotope_number)
|
|
{
|
|
coef = kit->second.Get_x_ratio_uncertainty();
|
|
|
|
/* scale epsilon in optimization equation */
|
|
|
|
my_array[((size_t)column - (size_t)col_epsilon) * max_column_count +
|
|
(size_t)column] = SCALE_EPSILON / coef;
|
|
|
|
/* set upper limit of change in positive direction */
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0;
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef;
|
|
sprintf(token, "%d%s %s",
|
|
(int) kit->second.Get_isotope_number(),
|
|
kit->second.Get_elt_name().c_str(), "eps+");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* set lower limit of change in negative direction */
|
|
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
my_array[count_rows * max_column_count + (size_t)i] = -coef;
|
|
sprintf(token, "%d%s %s",
|
|
(int) kit->second.Get_isotope_number(),
|
|
kit->second.Get_elt_name().c_str(), "eps-");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* inequalities for isotopes in phases
|
|
*/
|
|
/* row_isotope_phase_epsilon = count_rows; */
|
|
phase_isotope_inequalities(inv_ptr);
|
|
if (get_input_errors() > 0)
|
|
{
|
|
error_msg("Stopping because of input errors.", STOP);
|
|
}
|
|
/*
|
|
* Set non-negativity constraints
|
|
*/
|
|
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (inv_ptr->phases[i].constraint == PRECIPITATE)
|
|
{
|
|
delta[(size_t)col_phases + (size_t)i] = -1.0;
|
|
}
|
|
else if (inv_ptr->phases[i].constraint == DISSOLVE)
|
|
{
|
|
delta[(size_t)col_phases + (size_t)i] = 1.0;
|
|
}
|
|
}
|
|
for (i = 0; i < (inv_ptr->count_solns - 1); i++)
|
|
{
|
|
delta[i] = 1.0;
|
|
}
|
|
/*
|
|
* Scale water equation
|
|
*/
|
|
for (i = 0; i < max_column_count; i++)
|
|
{
|
|
my_array[(size_t)row_water * max_column_count + (size_t)i] *= SCALE_WATER;
|
|
}
|
|
/*
|
|
* Arrays are complete
|
|
*/
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
for (i = 0; i < count_unknowns; i++)
|
|
{
|
|
output_msg(sformatf( "%d\t%s\n", i, col_name[i]));
|
|
}
|
|
for (i = 0; i < count_rows; i++)
|
|
{
|
|
k = 0;
|
|
output_msg(sformatf( "%d\t%s\n", i, row_name[i]));
|
|
for (j = 0; j < count_unknowns + 1; j++)
|
|
{
|
|
if (k > 7)
|
|
{
|
|
output_msg(sformatf( "\n"));
|
|
k = 0;
|
|
}
|
|
output_msg(sformatf("%11.2e",
|
|
(double)my_array[(size_t)i * max_column_count + (size_t)j]));
|
|
k++;
|
|
}
|
|
if (k != 0)
|
|
{
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
output_msg(sformatf( "row_mb %d\n", row_mb));
|
|
output_msg(sformatf( "row_fract %d\n", row_fract));
|
|
output_msg(sformatf( "row_charge %d\n", row_charge));
|
|
output_msg(sformatf( "row_carbon %d\n", row_carbon));
|
|
output_msg(sformatf( "row_isotopes %d\n", row_isotopes));
|
|
output_msg(sformatf( "row_epsilon %d\n", row_epsilon));
|
|
|
|
output_msg(sformatf( "col_phases %d\n", col_phases));
|
|
output_msg(sformatf( "col_redox %d\n", col_redox));
|
|
output_msg(sformatf( "col_epsilon %d\n", col_epsilon));
|
|
output_msg(sformatf( "col_ph %d\n", col_ph));
|
|
output_msg(sformatf( "col_water %d\n", col_water));
|
|
output_msg(sformatf( "col_isotopes %d\n", col_isotopes));
|
|
output_msg(sformatf( "col_phase_isotopes %d\n",
|
|
col_phase_isotopes));
|
|
output_msg(sformatf( "count_unknowns %d\n", count_unknowns));
|
|
}
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
solve_inverse(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Exhaustively search for mass-balance models with two options
|
|
* -minimal on or off
|
|
* -range on or off
|
|
*
|
|
*/
|
|
int i, j, n;
|
|
int quit, print, first;
|
|
int first_of_model_size, model_size;
|
|
unsigned long minimal_bits, good_bits;
|
|
char token[MAX_LENGTH];
|
|
|
|
n = (int)count_unknowns; /* columns in A, C, E */
|
|
klmd = (max_row_count - 2);
|
|
nklmd = (n + klmd);
|
|
n2d = (size_t)n + 2;
|
|
|
|
max_good = MAX_MODELS;
|
|
max_bad = MAX_MODELS;
|
|
max_minimal = MAX_MODELS;
|
|
|
|
good.resize(max_good);
|
|
count_good = 0;
|
|
|
|
bad.resize(max_bad);
|
|
count_bad = 0;
|
|
|
|
minimal.resize(max_minimal);
|
|
count_minimal = 0;
|
|
|
|
col_back.resize(max_column_count);
|
|
row_back.resize(max_row_count);
|
|
|
|
/*
|
|
* Allocate space for arrays
|
|
*/
|
|
inv_cu.resize( 2 * (size_t)nklmd);
|
|
memset(&inv_cu[0], 0, ((2 * nklmd * sizeof(LDBLE))));
|
|
inv_iu.resize(2 * nklmd);
|
|
inv_is.resize(klmd);
|
|
|
|
for (i = 0; i < 79; i++)
|
|
token[i] = '=';
|
|
token[79] = '\0';
|
|
/*
|
|
* Set solutions, largest bit is final solution, smallest bit is initial solution 1
|
|
* Set phases, largest bit is last phase, smallest bit is first phase
|
|
* Set current bits to complete list.
|
|
*/
|
|
soln_bits = 0;
|
|
if (inv_ptr->count_solns + inv_ptr->phases.size() > 32)
|
|
{
|
|
error_msg
|
|
("For inverse modeling, sum of initial solutions and phases must be <= 32.\n\tFor all reasonable calculations, the sum should be much less than 32.",
|
|
STOP);
|
|
}
|
|
for (size_t i = inv_ptr->count_solns; i > 0; i--)
|
|
{
|
|
temp_bits = 1 << (i - 1);
|
|
soln_bits += temp_bits;
|
|
}
|
|
if (check_solns(inv_ptr) == ERROR)
|
|
{
|
|
error_msg("Calculations terminating.", STOP);
|
|
}
|
|
/*
|
|
* solutions are in highest bits, phases are in lower bits;
|
|
*/
|
|
/*
|
|
* All combinations of solutions
|
|
*/
|
|
first = TRUE;
|
|
for (; get_bits(soln_bits, (int)(inv_ptr->count_solns - 2),
|
|
(int)(inv_ptr->count_solns - 1)) > 0; soln_bits--)
|
|
{
|
|
/*
|
|
* Loop through all models of of descending size
|
|
*/
|
|
for (model_size = (int)inv_ptr->phases.size(); model_size >= 0; model_size--)
|
|
{
|
|
first_of_model_size = TRUE;
|
|
quit = TRUE;
|
|
while (next_set_phases(inv_ptr, first_of_model_size, model_size) == TRUE)
|
|
{
|
|
first_of_model_size = FALSE;
|
|
current_bits = (soln_bits << inv_ptr->phases.size()) + phase_bits;
|
|
|
|
if (subset_bad(current_bits) == TRUE
|
|
|| subset_minimal(current_bits) == TRUE)
|
|
continue;
|
|
quit = FALSE;
|
|
/*
|
|
* Switch for finding minimal models only
|
|
*/
|
|
if (inv_ptr->minimal == TRUE
|
|
&& superset_minimal(current_bits) == TRUE)
|
|
continue;
|
|
/*
|
|
* Solve for minimum epsilons, continue if no solution found.
|
|
*/
|
|
if (solve_with_mask(inv_ptr, current_bits) == ERROR)
|
|
{
|
|
save_bad(current_bits);
|
|
if (first == TRUE)
|
|
{
|
|
post_mortem();
|
|
quit = TRUE;
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
continue;
|
|
}
|
|
}
|
|
first = FALSE;
|
|
/*
|
|
* Model has been found, set bits
|
|
*/
|
|
good_bits = current_bits;
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (equal(inv_delta1[i + inv_ptr->count_solns], 0.0, TOL) == TRUE)
|
|
{
|
|
good_bits = set_bit(good_bits, (int)i, 0);
|
|
}
|
|
}
|
|
for (size_t i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, TOL) == TRUE)
|
|
{
|
|
good_bits = set_bit(good_bits, (int)(i + inv_ptr->phases.size()), 0);
|
|
}
|
|
}
|
|
/*
|
|
* Determine if model is new
|
|
*/
|
|
for (j = 0; j < count_good; j++)
|
|
{
|
|
if (good_bits == good[j])
|
|
break;
|
|
}
|
|
/*
|
|
* Calculate ranges and print model only if NOT looking for minimal models
|
|
*/
|
|
print = FALSE;
|
|
if (j >= count_good && inv_ptr->minimal == FALSE)
|
|
{
|
|
print = TRUE;
|
|
save_good(good_bits);
|
|
if (inv_ptr->range == TRUE)
|
|
{
|
|
range(inv_ptr, good_bits);
|
|
}
|
|
print_model(inv_ptr);
|
|
punch_model(inv_ptr);
|
|
dump_netpath_pat(inv_ptr);
|
|
}
|
|
/*
|
|
* If superset of a minimal model continue
|
|
*/
|
|
minimal_bits = good_bits;
|
|
if (superset_minimal(minimal_bits) == TRUE)
|
|
{
|
|
if (print == TRUE)
|
|
{
|
|
if (pr.inverse == TRUE && pr.all == TRUE)
|
|
{
|
|
output_msg(sformatf( "%s\n\n", token));
|
|
}
|
|
}
|
|
continue;
|
|
}
|
|
/*
|
|
* If not superset of minimal model, find minimal model
|
|
*/
|
|
minimal_bits = minimal_solve(inv_ptr, minimal_bits);
|
|
if (minimal_bits == good_bits && print == TRUE)
|
|
{
|
|
if (pr.inverse == TRUE && pr.all == TRUE)
|
|
{
|
|
output_msg(sformatf(
|
|
"\nModel contains minimum number of phases.\n"));
|
|
}
|
|
}
|
|
if (print == TRUE)
|
|
{
|
|
if (pr.inverse == TRUE && pr.all == TRUE)
|
|
{
|
|
output_msg(sformatf( "%s\n\n", token));
|
|
}
|
|
}
|
|
for (j = 0; j < count_good; j++)
|
|
{
|
|
if (minimal_bits == good[j])
|
|
break;
|
|
}
|
|
if (j >= count_good)
|
|
{
|
|
save_good(minimal_bits);
|
|
if (inv_ptr->range == TRUE)
|
|
{
|
|
range(inv_ptr, minimal_bits);
|
|
}
|
|
print_model(inv_ptr);
|
|
if (pr.inverse == TRUE && pr.all == TRUE)
|
|
{
|
|
output_msg(sformatf(
|
|
"\nModel contains minimum number of phases.\n"));
|
|
output_msg(sformatf( "%s\n\n", token));
|
|
}
|
|
punch_model(inv_ptr);
|
|
dump_netpath_pat(inv_ptr);
|
|
}
|
|
save_minimal(minimal_bits);
|
|
}
|
|
if (quit == TRUE)
|
|
break;
|
|
}
|
|
}
|
|
/*
|
|
* Summary print
|
|
*/
|
|
if (pr.inverse == TRUE && pr.all == TRUE)
|
|
{
|
|
output_msg(sformatf( "\nSummary of inverse modeling:\n\n"));
|
|
output_msg(sformatf( "\tNumber of models found: %d\n",
|
|
count_good));
|
|
output_msg(sformatf( "\tNumber of minimal models found: %d\n",
|
|
count_minimal));
|
|
output_msg(sformatf(
|
|
"\tNumber of infeasible sets of phases saved: %d\n",
|
|
count_bad));
|
|
output_msg(sformatf( "\tNumber of calls to cl1: %d\n",
|
|
count_calls));
|
|
}
|
|
my_array.clear();
|
|
delta.clear();
|
|
array1.clear();
|
|
inv_zero.clear();
|
|
inv_res.clear();
|
|
inv_delta1.clear();
|
|
delta2.clear();
|
|
delta3.clear();
|
|
delta_save.clear();
|
|
inv_cu.clear();
|
|
inv_iu.clear();
|
|
inv_is.clear();
|
|
col_name.clear();
|
|
row_name.clear();
|
|
col_back.clear();
|
|
row_back.clear();
|
|
min_delta.clear();
|
|
max_delta.clear();
|
|
good.clear();
|
|
bad.clear();
|
|
minimal.clear();
|
|
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
unsigned long Phreeqc::
|
|
minimal_solve(class inverse *inv_ptr, unsigned long minimal_bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Starting with phases indicated in minimal bits, sequentially
|
|
* remove phases to find minimal solution
|
|
*/
|
|
unsigned long temp_bits_l;
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "Beginning minimal solve: \n"));
|
|
bit_print(minimal_bits, (int)(inv_ptr->phases.size() + inv_ptr->count_solns));
|
|
}
|
|
for (size_t i = 0; i < inv_ptr->phases.size() + inv_ptr->count_solns - 1; i++)
|
|
{
|
|
if (get_bits(minimal_bits, (int)i, 1) == 0)
|
|
continue;
|
|
temp_bits_l = 1 << (int)i; /* 0's and one 1 */
|
|
temp_bits_l = ~temp_bits_l; /* 1's and one 0 */
|
|
minimal_bits = minimal_bits & temp_bits_l;
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "Solving for minimal\n"));
|
|
bit_print(minimal_bits, (int)(inv_ptr->phases.size() + inv_ptr->count_solns));
|
|
}
|
|
|
|
/*
|
|
* minimal_bits cannot be superset of a minimal model, but
|
|
* could be subset of one of the sets of minerals with no feasible solution
|
|
* If it is a subset, then replace mineral and go on to next
|
|
*/
|
|
if (subset_bad(minimal_bits) == TRUE)
|
|
{
|
|
/* put bit back */
|
|
minimal_bits = minimal_bits | ~temp_bits_l; /* 0's and one 1 */
|
|
continue;
|
|
}
|
|
if (solve_with_mask(inv_ptr, minimal_bits) == ERROR)
|
|
{
|
|
save_bad(minimal_bits);
|
|
/* put bit back */
|
|
minimal_bits = minimal_bits | ~temp_bits_l; /* 0's and one 1 */
|
|
}
|
|
|
|
}
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "\n\nMINIMAL MODEL\n\n"));
|
|
bit_print(minimal_bits, (int)(inv_ptr->phases.size() + inv_ptr->count_solns));
|
|
}
|
|
|
|
solve_with_mask(inv_ptr, minimal_bits);
|
|
unsigned long actual_bits = 0;
|
|
for (size_t i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, TOL) == FALSE)
|
|
{
|
|
actual_bits = set_bit(actual_bits, (int)(i + inv_ptr->phases.size()), 1);
|
|
}
|
|
}
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (equal(inv_delta1[i + inv_ptr->count_solns], 0.0, TOL) == FALSE)
|
|
{
|
|
actual_bits = set_bit(actual_bits, (int)i, 1);
|
|
}
|
|
}
|
|
if (actual_bits != minimal_bits)
|
|
{
|
|
warning_msg("Roundoff errors in minimal calculation");
|
|
}
|
|
return (actual_bits);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
solve_with_mask(class inverse *inv_ptr, unsigned long cur_bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Uses cur_bits to zero out columns of the array and then solves.
|
|
*/
|
|
int i, k, l, m, n;
|
|
|
|
/*
|
|
* Calculate dimensions
|
|
*/
|
|
k = (int)row_mb; /* rows in A */
|
|
l = (int)(row_epsilon - row_mb); /* rows in C */
|
|
m = (int)(count_rows - row_epsilon); /* rows in E */
|
|
n = (int)count_unknowns;
|
|
|
|
|
|
|
|
memcpy((void *) &(inv_res[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_row_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta2[0]), (void *) &(delta[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta_save[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
|
|
shrink(inv_ptr, &my_array[0], &array1[0],
|
|
&k, &l, &m, &n, cur_bits, &delta2[0], &col_back[0], &row_back[0]);
|
|
/*
|
|
* Save delta constraints
|
|
*/
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
delta_save[col_back[i]] = delta2[i];
|
|
}
|
|
|
|
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "\nColumns\n"));
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
output_msg(sformatf( "\t%d\t%s\n", i,
|
|
col_name[col_back[i]]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nRows\n"));
|
|
for (i = 0; i < k + l + m; i++)
|
|
{
|
|
output_msg(sformatf( "\t%d\t%s\n", i,
|
|
row_name[row_back[i]]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nA and B arrays:\n\n"));
|
|
array_print(&array1[0], k + l + m, n + 1, (int)max_column_count);
|
|
|
|
output_msg(sformatf( "\nInput delta vector:\n"));
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", i,
|
|
col_name[col_back[i]], (double) delta2[i]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
|
|
for (i = 0; i < k + l + m; i++)
|
|
{
|
|
if (inv_res[i] == 0)
|
|
continue;
|
|
output_msg(sformatf( "\nInput inv_res is non zero:\n"));
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", i,
|
|
row_name[row_back[i]], (double) inv_res[i]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
}
|
|
/*
|
|
* Call CL1
|
|
*/
|
|
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf(
|
|
"k, l, m, n, max_col, max_row\t%d\t%d\t%d\t%d\t%d\t%d\n",
|
|
k, l, m, n, max_column_count, max_row_count));
|
|
}
|
|
|
|
kode = 1;
|
|
iter = 100000;
|
|
count_calls++;
|
|
#ifdef INVERSE_CL1MP
|
|
if (inv_ptr->mp == TRUE)
|
|
{
|
|
cl1mp(k, l, m, n, (int)nklmd, (int)n2d, &array1[0],
|
|
&kode, inv_ptr->mp_tolerance, &iter, &delta2[0], &inv_res[0],
|
|
&error, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE, inv_ptr->mp_censor);
|
|
}
|
|
else
|
|
{
|
|
cl1(k, l, m, n, (int)nklmd, (int)n2d, &array1[0],
|
|
&kode, toler, &iter, &delta2[0], &inv_res[0],
|
|
&error, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE);
|
|
}
|
|
#else
|
|
cl1(k, l, m, n, (int)nklmd, (int)n2d, &array1[0],
|
|
&kode, toler, &iter, &delta2[0], &inv_res[0],
|
|
&error, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE);
|
|
#endif
|
|
if (kode == 3)
|
|
{
|
|
error_string = sformatf(
|
|
"Exceeded maximum iterations in inverse modeling: %d.\n"
|
|
"Recompile program with larger limit.", iter);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
memcpy((void *) &(inv_delta1[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
inv_delta1[col_back[i]] = delta2[i];
|
|
}
|
|
|
|
/*
|
|
* Debug, write results
|
|
*/
|
|
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "kode: %d\titer: %d\terror: %e\n", kode,
|
|
iter, (double) error));
|
|
output_msg(sformatf( "\nsolution vector:\n"));
|
|
for (i = 0; i < n; i++)
|
|
{
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", i,
|
|
col_name[col_back[i]], (double) delta2[i]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
|
|
output_msg(sformatf( "\nresidual vector:\n"));
|
|
for (i = 0; i < (k + l + m); i++)
|
|
{
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e\n", i,
|
|
row_name[row_back[i]], (double) inv_res[i]));
|
|
}
|
|
}
|
|
|
|
if (kode != 0)
|
|
{
|
|
return (ERROR);
|
|
}
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
unsigned long Phreeqc::
|
|
get_bits(unsigned long bits, int position, int number)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Returns number of bits from position and below.
|
|
* position begins at 0.
|
|
*/
|
|
return ((bits >> (position + 1 - number)) & ~(~0ul << number));
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
save_minimal(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Keeps list of minimal models
|
|
*/
|
|
minimal[count_minimal] = bits;
|
|
count_minimal++;
|
|
if (count_minimal >= max_minimal)
|
|
{
|
|
max_minimal *= 2;
|
|
minimal.resize(max_minimal);
|
|
}
|
|
return (TRUE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
save_good(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Keeps list of good models, not necessarily minimal
|
|
*/
|
|
good[count_good] = bits;
|
|
count_good++;
|
|
if (count_good >= max_good)
|
|
{
|
|
max_good *= 2;
|
|
good.resize(max_good);
|
|
}
|
|
return (TRUE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
save_bad(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Keeps list of sets of phases with no feasible solution
|
|
*/
|
|
bad[count_bad] = bits;
|
|
count_bad++;
|
|
if (count_bad >= max_bad)
|
|
{
|
|
max_bad *= 2;
|
|
bad.resize(max_bad);
|
|
}
|
|
return (TRUE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
superset_minimal(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Checks whether bits is a superset of any of the minimal models
|
|
*/
|
|
int i;
|
|
unsigned long temp_bits_l;
|
|
for (i = 0; i < count_minimal; i++)
|
|
{
|
|
temp_bits_l = bits | minimal[i];
|
|
if (temp_bits_l == bits)
|
|
{
|
|
return (TRUE);
|
|
}
|
|
}
|
|
return (FALSE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
subset_bad(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Checks whether bits is a superset of any of the bad models
|
|
*/
|
|
int i;
|
|
unsigned long temp_bits_l;
|
|
for (i = 0; i < count_bad; i++)
|
|
{
|
|
temp_bits_l = bits | bad[i];
|
|
if (temp_bits_l == bad[i])
|
|
{
|
|
return (TRUE);
|
|
}
|
|
}
|
|
return (FALSE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
subset_minimal(unsigned long bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Checks whether bits is a subset of any of the minimal models
|
|
*/
|
|
int i;
|
|
unsigned long temp_bits_l;
|
|
for (i = 0; i < count_minimal; i++)
|
|
{
|
|
temp_bits_l = bits | minimal[i];
|
|
if (temp_bits_l == minimal[i])
|
|
{
|
|
return (TRUE);
|
|
}
|
|
}
|
|
return (FALSE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
bit_print(unsigned long bits, int l)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Prints l bits of an unsigned long
|
|
*/
|
|
int i;
|
|
|
|
for (i = l - 1; i >= 0; i--)
|
|
{
|
|
output_msg(sformatf( "%lu ", get_bits(bits, i, 1)));
|
|
}
|
|
output_msg(sformatf( "\n"));
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
print_model(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Prints model
|
|
*/
|
|
int i, j;
|
|
size_t column;
|
|
int print_msg;
|
|
cxxSolution *solution_ptr;
|
|
class master *master_ptr;
|
|
LDBLE d1, d2, d3, d4;
|
|
char token[MAX_LENGTH];
|
|
/*
|
|
* Update screen
|
|
*/
|
|
status(count_good, NULL);
|
|
/*
|
|
* print solution data, epsilons, and revised data
|
|
*/
|
|
if (pr.inverse == FALSE || pr.all == FALSE)
|
|
return (OK);
|
|
max_pct = 0;
|
|
scaled_error = 0;
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, toler) == TRUE)
|
|
continue;
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
xsolution_zero();
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
master_ptr = master_bsearch(jit->first.c_str());
|
|
master_ptr->total = jit->second;
|
|
}
|
|
|
|
output_msg(sformatf( "\nSolution %d: %s\n", inv_ptr->solns[i],
|
|
solution_ptr->Get_description().c_str()));
|
|
output_msg(sformatf(
|
|
"\n%15.15s %12.12s %12.12s %12.12s\n", " ",
|
|
"Input", "Delta", "Input+Delta"));
|
|
master_alk->total = solution_ptr->Get_total_alkalinity();
|
|
if (inv_ptr->carbon == TRUE)
|
|
{
|
|
d1 = solution_ptr->Get_ph();
|
|
d2 = inv_delta1[col_ph + i] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
output_msg(sformatf(
|
|
"%15.15s %12.3e +%12.3e =%12.3e\n", "pH",
|
|
(double) d1, (double) d2, (double) d3));
|
|
if (inv_ptr->ph_uncertainties[i] > 0)
|
|
{
|
|
scaled_error += fabs(d2) / inv_ptr->ph_uncertainties[i];
|
|
/* debug
|
|
output_msg(sformatf( "%e\t%e\t%e\n", fabs(d2) / inv_ptr->ph_uncertainties[i], fabs(d2), inv_ptr->ph_uncertainties[i]));
|
|
*/
|
|
}
|
|
else if (d2 != 0.0)
|
|
{
|
|
error_msg("Computing delta pH/uncertainty", CONTINUE);
|
|
}
|
|
}
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
if (inv_ptr->elts[j].master->s == s_eminus)
|
|
continue;
|
|
d1 = inv_ptr->elts[j].master->total;
|
|
d2 = inv_delta1[col_epsilon + j * inv_ptr->count_solns +
|
|
i] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
|
|
output_msg(sformatf(
|
|
"%15.15s %12.3e +%12.3e =%12.3e\n",
|
|
inv_ptr->elts[j].master->elt->name, (double) d1,
|
|
(double) d2, (double) d3));
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == FALSE)
|
|
{
|
|
d3 = fabs(d2 / d1);
|
|
if (d3 > max_pct)
|
|
max_pct = d3;
|
|
}
|
|
d4 = 0;
|
|
if (inv_ptr->elts[j].uncertainties[i] > 0)
|
|
{
|
|
d4 = fabs(inv_ptr->elts[j].uncertainties[i] * d1);
|
|
}
|
|
else if (inv_ptr->elts[j].uncertainties[i] < 0)
|
|
{
|
|
d4 = -inv_ptr->elts[j].uncertainties[i];
|
|
}
|
|
if (d4 > 0)
|
|
{
|
|
scaled_error += fabs(d2) / d4;
|
|
/* debug
|
|
output_msg(sformatf( "%e\t%e\t%e\n", fabs(d2) / d4, fabs(d2), d4));
|
|
*/
|
|
}
|
|
else if (d2 != 0.0)
|
|
{
|
|
error_msg("Computing delta element/uncertainty", CONTINUE);
|
|
}
|
|
}
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
/* adjustments to solution isotope composition */
|
|
for (j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (inv_ptr->isotope_unknowns[j].elt_name !=
|
|
string_hsave(kit->second.Get_elt_name().c_str()) ||
|
|
inv_ptr->isotope_unknowns[j].isotope_number !=
|
|
kit->second.Get_isotope_number())
|
|
continue;
|
|
d1 = kit->second.Get_ratio();
|
|
d2 = inv_delta1[col_isotopes +
|
|
i * inv_ptr->isotope_unknowns.size() +
|
|
j] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
sprintf(token, "%d%s",
|
|
(int) inv_ptr->isotope_unknowns[j].
|
|
isotope_number,
|
|
inv_ptr->isotope_unknowns[j].elt_name);
|
|
output_msg(sformatf(
|
|
"%15.15s %12g +%12g =%12g\n", token,
|
|
(double) d1, (double) d2, (double) d3));
|
|
/*
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == FALSE ) {
|
|
d3 = fabs(d2/d1);
|
|
if (d3 > max_pct) max_pct = d3;
|
|
}
|
|
*/
|
|
if (kit->second.Get_x_ratio_uncertainty() > 0)
|
|
{
|
|
scaled_error +=
|
|
fabs(d2) /
|
|
kit->second.Get_x_ratio_uncertainty();
|
|
/* debug
|
|
output_msg(sformatf( "%e\t%e\t%e\n", fabs(d2) / solution_ptr->isotopes[k].x_ratio_uncertainty , fabs(d2), solution_ptr->isotopes[k].x_ratio_uncertainty));
|
|
*/
|
|
}
|
|
else if (d2 != 0.0)
|
|
{
|
|
error_msg
|
|
("Computing delta solution isotope/uncertainty",
|
|
CONTINUE);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Adjustments to phases
|
|
*/
|
|
print_msg = FALSE;
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
output_msg(sformatf( "\nIsotopic composition of phases:\n"));
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (inv_ptr->phases[i].isotopes.size() == 0)
|
|
continue;
|
|
size_t j = col_phases + i;
|
|
if (equal(inv_delta1[j], 0.0, toler) == TRUE &&
|
|
equal(min_delta[j], 0.0, toler) == TRUE &&
|
|
equal(max_delta[j], 0.0, toler) == TRUE)
|
|
continue;
|
|
std::vector<class isotope>& isotope_ref = inv_ptr->phases[i].isotopes;
|
|
for (size_t j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
for (size_t k = 0; k < inv_ptr->phases[i].isotopes.size(); k++)
|
|
{
|
|
if (inv_ptr->isotopes[j].elt_name !=
|
|
isotope_ref[k].elt_name ||
|
|
inv_ptr->isotopes[j].isotope_number !=
|
|
isotope_ref[k].isotope_number)
|
|
continue;
|
|
d1 = isotope_ref[k].ratio;
|
|
column = col_phase_isotopes + i * inv_ptr->isotopes.size() + j;
|
|
if (inv_delta1[col_phases + i] != 0.0)
|
|
{
|
|
d2 = inv_delta1[column] / inv_delta1[col_phases + i];
|
|
}
|
|
else
|
|
{
|
|
continue;
|
|
}
|
|
d3 = d1 + d2;
|
|
if (equal(d1, 0.0, 1e-7) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, 1e-7) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, 1e-7) == TRUE)
|
|
d3 = 0.0;
|
|
sprintf(token, "%d%s %s",
|
|
(int) inv_ptr->isotopes[j].isotope_number,
|
|
inv_ptr->isotopes[j].elt_name,
|
|
inv_ptr->phases[i].phase->name);
|
|
output_msg(sformatf(
|
|
"%15.15s %12g +%12g =%12g", token,
|
|
(double) d1, (double) d2, (double) d3));
|
|
if (fabs(d2) > (isotope_ref[k].ratio_uncertainty + toler))
|
|
{
|
|
output_msg(sformatf( " **"));
|
|
print_msg = TRUE;
|
|
}
|
|
output_msg(sformatf( "\n"));
|
|
if (isotope_ref[k].ratio_uncertainty > 0)
|
|
{
|
|
scaled_error +=
|
|
fabs(d2) / isotope_ref[k].ratio_uncertainty;
|
|
/* debug
|
|
output_msg(sformatf( "%e\t%e\t%e\n", fabs(d2) / isotope_ptr[k].ratio_uncertainty, fabs(d2), isotope_ptr[k].ratio_uncertainty));
|
|
*/
|
|
}
|
|
else if (d2 != 0.0)
|
|
{
|
|
error_msg
|
|
("Computing delta phase isotope/uncertainty",
|
|
CONTINUE);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (print_msg == TRUE)
|
|
{
|
|
output_msg(sformatf(
|
|
"\n**\tWARNING: The adjustment to at least one isotopic"
|
|
"\n\tcomposition of a phase exceeded the specified uncertainty"
|
|
"\n\tfor the phase. If the phase is not constrained to dissolve"
|
|
"\n\tor precipitate, then the isotopic composition of the phase"
|
|
"\n\tis also unconstrained.\n"));
|
|
}
|
|
output_msg(sformatf( "\n%-20.20s %7s %12.12s %12.12s\n",
|
|
"Solution fractions:", " ", "Minimum", "Maximum"));
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
d1 = inv_delta1[i];
|
|
d2 = min_delta[i];
|
|
d3 = max_delta[i];
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
output_msg(sformatf( "%11s%4d %12.3e %12.3e %12.3e\n",
|
|
"Solution", inv_ptr->solns[i], (double) d1, (double) d2,
|
|
(double) d3));
|
|
}
|
|
|
|
// appt, calculate and print SI's
|
|
LDBLE t_i, p_i, iap, lk, t;
|
|
const char *name;
|
|
class rxn_token *rxn_ptr;
|
|
CReaction *reaction_ptr;
|
|
|
|
output_msg(sformatf( "\n%-25.25s %2s %12.12s %12.12s %-18.18s (Approximate SI in solution ",
|
|
"Phase mole transfers:", " ", "Minimum", "Maximum", "Formula"));
|
|
|
|
for (i = 0; i < inv_ptr->count_solns - 1; i++)
|
|
output_msg(sformatf("%d, ", inv_ptr->solns[i]));
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
t_i = solution_ptr->Get_tc() + 273.15;
|
|
p_i = solution_ptr->Get_patm();
|
|
output_msg(sformatf("%d at %3d K, %3d atm)\n", inv_ptr->solns[i], int(t_i), int(floor(p_i + 0.5))));
|
|
p_i *= PASCAL_PER_ATM;
|
|
|
|
for (size_t i = col_phases; i < col_redox; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, toler) == TRUE &&
|
|
equal(min_delta[i], 0.0, toler) == TRUE &&
|
|
equal(max_delta[i], 0.0, toler) == TRUE)
|
|
continue;
|
|
|
|
d1 = inv_delta1[i];
|
|
d2 = min_delta[i];
|
|
d3 = max_delta[i];
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
output_msg(sformatf(
|
|
"%15.15s %12.3e %12.3e %12.3e %-25.25s (", col_name[i],
|
|
(double)d1, (double)d2, (double)d3, inv_ptr->phases[i - col_phases].phase->formula));
|
|
|
|
size_t i1 = 0;
|
|
for (; i1 < phases.size(); i1++)
|
|
{
|
|
if (Utilities::strcmp_nocase(phases[i1]->name, col_name[i]))
|
|
continue;
|
|
reaction_ptr = &phases[i1]->rxn_s;
|
|
for (size_t i2 = 0; i2 < inv_ptr->count_solns; i2++)
|
|
{
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i2]);
|
|
|
|
reaction_ptr->logk[delta_v] = calc_delta_v(*reaction_ptr, true) - phases[i1]->logk[vm0];
|
|
if (reaction_ptr->logk[delta_v])
|
|
mu_terms_in_logk = true;
|
|
lk = k_calc(reaction_ptr->logk, t_i, p_i);
|
|
|
|
iap = 0.0;
|
|
for (rxn_ptr = &reaction_ptr->token[0] + 1; rxn_ptr->s != NULL; rxn_ptr++)
|
|
{
|
|
t = 0;
|
|
if (rxn_ptr->s == s_eminus)
|
|
t = -solution_ptr->Get_pe();
|
|
else if (!Utilities::strcmp_nocase(rxn_ptr->s->name, "H2O"))
|
|
t = log10(solution_ptr->Get_ah2o());
|
|
else if (!Utilities::strcmp_nocase(rxn_ptr->s->name, "H+"))
|
|
t = -solution_ptr->Get_ph();
|
|
else
|
|
{
|
|
if (rxn_ptr->s->secondary)
|
|
name = rxn_ptr->s->secondary->elt->name;
|
|
else
|
|
name = rxn_ptr->s->primary->elt->name;
|
|
t = solution_ptr->Get_master_activity()[name];
|
|
}
|
|
if (t)
|
|
iap += t * rxn_ptr->coef;
|
|
else
|
|
{
|
|
iap = -999; break;
|
|
}
|
|
}
|
|
|
|
if (iap == -999)
|
|
output_msg(sformatf(" "));
|
|
else
|
|
output_msg(sformatf("%6.2f", iap - lk));
|
|
if (i2 < inv_ptr->count_solns - 1)
|
|
output_msg(sformatf(","));
|
|
}
|
|
}
|
|
output_msg(sformatf(")\n"));
|
|
}
|
|
output_msg(sformatf( "\n%-25.25s\n", "Redox mole transfers:"));
|
|
for (size_t i = col_redox; i < col_epsilon; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, toler) == TRUE)
|
|
continue;
|
|
output_msg(sformatf( "%15.15s %12.3e\n", col_name[i],
|
|
(double) inv_delta1[i]));
|
|
}
|
|
|
|
output_msg(sformatf(
|
|
"\nSum of residuals (epsilons in documentation): %12.3e\n",
|
|
((double) (error / SCALE_EPSILON))));
|
|
output_msg(sformatf(
|
|
"Sum of delta/uncertainty limit: %12.3e\n",
|
|
(double) scaled_error));
|
|
output_msg(sformatf(
|
|
"Maximum fractional error in element concentration: %12.3e\n",
|
|
(double) max_pct));
|
|
/*
|
|
* Flush buffer after each model
|
|
*/
|
|
output_flush();
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
punch_model_heading(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Prints model headings to selected output file
|
|
*/
|
|
int i;
|
|
char token[MAX_LENGTH];
|
|
|
|
std::map < int, SelectedOutput >::iterator so_it = SelectedOutput_map.begin();
|
|
for ( ; so_it != SelectedOutput_map.end(); so_it++)
|
|
{
|
|
// set punch file
|
|
current_selected_output = &(so_it->second);
|
|
if (pr.punch == FALSE ||
|
|
current_selected_output == NULL ||
|
|
!current_selected_output->Get_inverse() ||
|
|
!current_selected_output->Get_active())
|
|
continue;
|
|
phrq_io->Set_punch_ostream(current_selected_output->Get_punch_ostream());
|
|
|
|
int l = (!current_selected_output->Get_high_precision()) ? 15 : 20;
|
|
inverse_heading_names.clear();
|
|
/*
|
|
* Print sum of residuals and maximum fractional error
|
|
*/
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, "Sum_resid"));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, "Sum_Delta/U"));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, "MaxFracErr"));
|
|
|
|
/*
|
|
* Print solution numbers
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
sprintf(token, "Soln_%d", inv_ptr->solns[i]);
|
|
std::string tok1(token);
|
|
tok1.append("_min");
|
|
std::string tok2(token);
|
|
tok2.append("_max");
|
|
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, token));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, tok1.c_str()));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, tok2.c_str()));
|
|
}
|
|
/*
|
|
* Print phase names
|
|
*/
|
|
for (size_t i = col_phases; i < col_redox; i++)
|
|
{
|
|
|
|
std::string tok1(col_name[i]);
|
|
tok1.append("_min");
|
|
std::string tok2(col_name[i]);
|
|
tok2.append("_max");
|
|
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, col_name[i]));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, tok1.c_str()));
|
|
inverse_heading_names.push_back(sformatf("%*s\t", l, tok2.c_str()));
|
|
}
|
|
for (size_t j = 0; j < inverse_heading_names.size(); j++)
|
|
{
|
|
fpunchf_heading(inverse_heading_names[j].c_str());
|
|
//user_punch_headings[j] = string_hsave(heading_names[j].c_str());
|
|
}
|
|
fpunchf_heading("\n");
|
|
}
|
|
current_selected_output = NULL;
|
|
phrq_io->Set_punch_ostream(NULL);
|
|
/*
|
|
* Flush buffer after each model
|
|
*/
|
|
punch_flush();
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
punch_model(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Prints model to selected output file
|
|
*/
|
|
int i;
|
|
LDBLE d1, d2, d3;
|
|
//if (punch.in == FALSE || pr.punch == FALSE || punch.inverse == FALSE)
|
|
// return (OK);
|
|
|
|
UserPunch temp_user_punch;
|
|
current_user_punch = & temp_user_punch;
|
|
temp_user_punch.Set_headings(inverse_heading_names);
|
|
|
|
std::map < int, SelectedOutput >::iterator so_it = SelectedOutput_map.begin();
|
|
for ( ; so_it != SelectedOutput_map.end(); so_it++)
|
|
{
|
|
current_selected_output = &(so_it->second);
|
|
if (pr.punch == FALSE ||
|
|
current_selected_output == NULL ||
|
|
!current_selected_output->Get_inverse() ||
|
|
!current_selected_output->Get_active())
|
|
continue;
|
|
phrq_io->Set_punch_ostream(current_selected_output->Get_punch_ostream());
|
|
|
|
n_user_punch_index = 0;
|
|
/*
|
|
* write residual info
|
|
*/
|
|
if (!current_selected_output->Get_high_precision())
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) (error / SCALE_EPSILON));
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) scaled_error);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) max_pct);
|
|
}
|
|
else
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) (error / SCALE_EPSILON));
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) scaled_error);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) max_pct);
|
|
}
|
|
/*
|
|
* write solution fractions
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
d1 = inv_delta1[i];
|
|
d2 = min_delta[i];
|
|
d3 = max_delta[i];
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
if (!current_selected_output->Get_high_precision())
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d1);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d2);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d3);
|
|
}
|
|
else
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d1);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d2);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d3);
|
|
}
|
|
}
|
|
/*
|
|
* write phase transfers
|
|
*/
|
|
for (size_t i = col_phases; i < col_redox; i++)
|
|
{
|
|
d1 = inv_delta1[i];
|
|
d2 = min_delta[i];
|
|
d3 = max_delta[i];
|
|
if (equal(d1, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d1 = 0.0;
|
|
if (equal(d2, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d2 = 0.0;
|
|
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
|
|
d3 = 0.0;
|
|
if (!current_selected_output->Get_high_precision())
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d1);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d2);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%15.4e\t", (double) d3);
|
|
}
|
|
else
|
|
{
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d1);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d2);
|
|
fpunchf(trim(inverse_heading_names[n_user_punch_index++]).c_str(), "%20.12e\t", (double) d3);
|
|
}
|
|
}
|
|
punch_msg("\n");
|
|
|
|
/*
|
|
* Flush buffer after each model
|
|
*/
|
|
punch_flush();
|
|
}
|
|
current_selected_output = NULL;
|
|
phrq_io->Set_punch_ostream(NULL);
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
unsigned long Phreeqc::
|
|
set_bit(unsigned long bits, int position, int value)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Sets a single bit
|
|
*/
|
|
unsigned long temp_bits_l;
|
|
|
|
temp_bits_l = 1 << position;
|
|
if (value == 0)
|
|
{
|
|
temp_bits_l = ~temp_bits_l;
|
|
temp_bits_l = bits & temp_bits_l;
|
|
}
|
|
else
|
|
{
|
|
temp_bits_l = bits | temp_bits_l;
|
|
}
|
|
return (temp_bits_l);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
next_set_phases(class inverse *inv_ptr,
|
|
int first_of_model_size, int model_size)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
int i, j, k;
|
|
unsigned long temp_bits_l;
|
|
|
|
/*
|
|
* min_ and max_position are arrays, logically with length
|
|
* of model_size, that contain minimum and maximum
|
|
* phase numbers that can be in that model position.
|
|
*
|
|
* now contains a list of phase numbers to mark as in for this
|
|
* model
|
|
*/
|
|
|
|
/*
|
|
* Initialize for a given model_size
|
|
*/
|
|
if (first_of_model_size == TRUE)
|
|
{
|
|
for (i = 0; i < model_size; i++)
|
|
{
|
|
min_position[i] = i;
|
|
now[i] = i;
|
|
max_position[i] = (int)inv_ptr->phases.size() - model_size + i;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
/*
|
|
* Determine next combination of phases for fixed model_size
|
|
*/
|
|
for (i = (model_size - 1); i >= 0; i--)
|
|
{
|
|
if (now[i] < max_position[i])
|
|
{
|
|
now[i]++;
|
|
if (i < (model_size - 1))
|
|
{
|
|
k = now[i];
|
|
for (j = (i + 1); j < model_size; j++)
|
|
{
|
|
k++;
|
|
now[j] = k;
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
if (i < 0)
|
|
return (FALSE);
|
|
}
|
|
/*
|
|
* Set bits which switch in phases
|
|
*/
|
|
temp_bits_l = 0;
|
|
for (j = 0; j < model_size; j++)
|
|
{
|
|
temp_bits_l += (1 << now[j]);
|
|
}
|
|
phase_bits = temp_bits_l;
|
|
return (TRUE);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
range(class inverse *inv_ptr, unsigned long cur_bits)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Takes the model from cur_bits and sequentially determines the
|
|
* minimum and maximum values for each solution fraction and
|
|
* each phase mass transfer.
|
|
*/
|
|
int i, j;
|
|
int k, l, m, n;
|
|
int f;
|
|
unsigned long bits;
|
|
LDBLE error2;
|
|
/*
|
|
* Include forced solutions and phases in range calculation
|
|
*/
|
|
for (size_t i = 0; i < inv_ptr->count_solns + inv_ptr->phases.size(); i++)
|
|
{
|
|
if (i < inv_ptr->phases.size())
|
|
{
|
|
if (inv_ptr->phases[i].force == TRUE)
|
|
{
|
|
cur_bits = set_bit(cur_bits, (int)i, 1);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (inv_ptr->force_solns[i - inv_ptr->phases.size()] == TRUE)
|
|
{
|
|
cur_bits = set_bit(cur_bits, (int)i, 1);
|
|
}
|
|
}
|
|
}
|
|
|
|
memcpy((void *) &(min_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(max_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
/*
|
|
* Switch bits so that phases are high and solutions are low
|
|
*/
|
|
bits =
|
|
get_bits(cur_bits, (int)(inv_ptr->phases.size() + inv_ptr->count_solns) - 1,
|
|
(int)inv_ptr->count_solns);
|
|
bits +=
|
|
(get_bits(cur_bits, (int)inv_ptr->phases.size() - 1, (int)inv_ptr->phases.size())
|
|
<< (int)inv_ptr->count_solns);
|
|
/*
|
|
* Do range calculation
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns + inv_ptr->phases.size(); i++)
|
|
{
|
|
if (inv_ptr->count_solns == i + 1)
|
|
{
|
|
min_delta[i] = 1.0;
|
|
max_delta[i] = 1.0;
|
|
continue;
|
|
}
|
|
if (get_bits(bits, i, 1) == 0)
|
|
continue;
|
|
/*
|
|
* Calculate min and max
|
|
*/
|
|
for (f = -1; f < 2; f += 2)
|
|
{
|
|
k = (int)row_mb; /* rows in A */
|
|
l = (int)(row_epsilon - row_mb); /* rows in C */
|
|
m = (int)(count_rows - row_epsilon); /* rows in E */
|
|
n = (int)count_unknowns; /* number of variables */
|
|
/*
|
|
* Copy equations
|
|
*/
|
|
memcpy((void *) &(array1[0]), (void *) &(my_array[0]),
|
|
(size_t) max_column_count * max_row_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta2[0]), (void *) &(delta[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta3[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta_save[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(inv_res[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_row_count * sizeof(LDBLE));
|
|
|
|
/*
|
|
* Change optimization
|
|
*/
|
|
for (j = 0; j < k; j++)
|
|
{
|
|
memcpy((void *) &(array1[j * max_column_count]),
|
|
(void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
}
|
|
array1[i] = 1.0;
|
|
if (f < 1)
|
|
{
|
|
array1[n] = -fabs(inv_ptr->range_max);
|
|
}
|
|
else
|
|
{
|
|
array1[n] = fabs(inv_ptr->range_max);
|
|
}
|
|
shrink(inv_ptr, &array1[0], &array1[0],
|
|
&k, &l, &m, &n, cur_bits, &delta2[0], &col_back[0], &row_back[0]);
|
|
/*
|
|
* Save delta constraints
|
|
*/
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
delta_save[col_back[j]] = delta2[j];
|
|
}
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "\nInput delta:\n\n"));
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
output_msg(sformatf( "\t%d %s\t%g\n", j,
|
|
col_name[col_back[j]], (double) delta2[j]));
|
|
}
|
|
output_msg(sformatf( "\nA and B arrays:\n\n"));
|
|
array_print(&array1[0], k + l + m, n + 1, (int)max_column_count);
|
|
}
|
|
kode = 1;
|
|
iter = 200;
|
|
count_calls++;
|
|
#ifdef INVERSE_CL1MP
|
|
if (inv_ptr->mp == TRUE)
|
|
{
|
|
cl1mp(k, l, m, n, (int)nklmd, (int)n2d, &array1[0],
|
|
&kode, inv_ptr->mp_tolerance, &iter, &delta2[0], &inv_res[0],
|
|
&error2, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE, inv_ptr->mp_censor);
|
|
}
|
|
else
|
|
{
|
|
cl1(k, l, m, n, (int)nklmd, (int)n2d, &array1[0], &kode, toler, &iter, &delta2[0],
|
|
&inv_res[0], &error2, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE);
|
|
}
|
|
#else
|
|
cl1(k, l, m, n, (int)nklmd, (int)n2d, &array1[0], &kode, toler, &iter, &delta2[0],
|
|
&inv_res[0], &error2, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE);
|
|
#endif
|
|
if (kode != 0)
|
|
{
|
|
output_msg(sformatf(
|
|
"Error in subroutine range. Kode = %d\n", kode));
|
|
}
|
|
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "kode: %d\titer: %d\terror: %e\n",
|
|
kode, iter, (double) error2));
|
|
output_msg(sformatf( "k, l, m, n: %d\t%d\t%d\t%d\n", k,
|
|
l, m, n));
|
|
output_msg(sformatf( "\nsolution vector %s\n",
|
|
col_name[i]));
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", j,
|
|
col_name[col_back[j]], (double) delta2[j]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
|
|
output_msg(sformatf( "\nresidual vector:\n"));
|
|
for (j = 0; j < (k + l + m); j++)
|
|
{
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e\n", j,
|
|
row_name[row_back[j]], (double) inv_res[j]));
|
|
}
|
|
}
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
if (col_back[j] == i)
|
|
break;
|
|
}
|
|
if (f < 0)
|
|
{
|
|
min_delta[i] = delta2[j];
|
|
}
|
|
else
|
|
{
|
|
max_delta[i] = delta2[j];
|
|
}
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
delta3[col_back[j]] = delta2[j];
|
|
}
|
|
}
|
|
}
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
shrink(class inverse *inv_ptr, LDBLE * array_in, LDBLE * array_out,
|
|
int *k, int *l, int *m, int *n,
|
|
unsigned long cur_bits,
|
|
LDBLE * delta_l, int *col_back_l, int *row_back_l)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Shrink eliminates any rows that are all zeros and any columns
|
|
* that are not in cur_bits, result is put in array_out
|
|
*
|
|
* k, l, m, n return the new sizes of the array.
|
|
* delta is remapped to retain any non-negativity constraints
|
|
* Col_back maps columns that remain back to original columns
|
|
* Row_back maps rows that remain back to original rows
|
|
*/
|
|
int i, j, row;
|
|
int k1, l1, m1;
|
|
size_t cur_col, column;
|
|
int nonzero;
|
|
/*
|
|
* Copy array_in to array_out
|
|
*/
|
|
if (array_in != array_out)
|
|
{
|
|
for (i = 0; i < (*k + *l + *m); i++)
|
|
{
|
|
memcpy(&(array_out[i * max_column_count]),
|
|
&(array_in[i * max_column_count]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
}
|
|
}
|
|
/*
|
|
* Determine columns to eliminate
|
|
*/
|
|
for (i = 0; i < (*n + 1); i++)
|
|
col_back_l[i] = i;
|
|
|
|
/*
|
|
* Drop phases not in model
|
|
*/
|
|
for (i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (get_bits(cur_bits, i, 1) == 0)
|
|
{
|
|
col_back_l[col_phases + i] = -1;
|
|
/* drop isotopes */
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
for (j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
column =
|
|
col_phase_isotopes + i * inv_ptr->isotopes.size() + j;
|
|
col_back_l[column] = -1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Drop solutions not in model
|
|
*/
|
|
for (i = 0; i < (inv_ptr->count_solns - 1); i++)
|
|
{
|
|
if (get_bits(cur_bits, (int)inv_ptr->phases.size() + i, 1) == 0)
|
|
{
|
|
col_back_l[i] = -1;
|
|
/* drop all epsilons for the solution */
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
column = col_epsilon + j * inv_ptr->count_solns + i;
|
|
col_back_l[column] = -1;
|
|
}
|
|
/* drop pH for the solution */
|
|
if (inv_ptr->carbon == TRUE)
|
|
{
|
|
column = col_ph + i;
|
|
col_back_l[column] = -1;
|
|
}
|
|
/* drop isotopes */
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
for (size_t j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
column = col_isotopes +
|
|
i * inv_ptr->isotope_unknowns.size() + j;
|
|
col_back_l[column] = -1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Drop epsilons not used
|
|
*/
|
|
for (i = (int)col_epsilon; i < *n; i++)
|
|
{
|
|
if (col_back_l[i] < 0)
|
|
continue;
|
|
for (j = 0; j < (*k + *l + *m); j++)
|
|
{
|
|
if (array_out[(size_t)j * max_column_count + (size_t)i] != 0)
|
|
break;
|
|
}
|
|
if (j == (*k + *l + *m))
|
|
{
|
|
col_back_l[i] = -1;
|
|
}
|
|
}
|
|
/*
|
|
* rewrite array_out
|
|
*/
|
|
cur_col = 0;
|
|
for (i = 0; i < (*n + 1); i++)
|
|
{
|
|
if (col_back_l[i] < 0)
|
|
continue;
|
|
if (cur_col == col_back_l[i])
|
|
{
|
|
cur_col++;
|
|
continue;
|
|
}
|
|
for (j = 0; j < (*k + *l + *m); j++)
|
|
{
|
|
array_out[j * max_column_count + cur_col] =
|
|
array_out[j * max_column_count + i];
|
|
}
|
|
col_back_l[cur_col] = col_back_l[i];
|
|
delta_l[cur_col] = delta_l[i];
|
|
cur_col++;
|
|
}
|
|
*n = (int)cur_col - 1;
|
|
/*
|
|
* Eliminate unnecessary optimization eqns
|
|
*/
|
|
row = 0;
|
|
k1 = 0;
|
|
for (i = 0; i < *k; i++)
|
|
{
|
|
if (memcmp(&(array_out[i * max_column_count]), &(inv_zero[0]),
|
|
(size_t) (*n) * sizeof(LDBLE)) == 0)
|
|
{
|
|
continue;
|
|
}
|
|
/*
|
|
memcpy(&(array_out[row * max_column_count]), &(array_out[i * max_column_count]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
*/
|
|
if (i > row)
|
|
{
|
|
if ((row*max_column_count + (*n + 1)) >= (i * max_column_count))
|
|
{
|
|
assert(false);
|
|
}
|
|
|
|
|
|
memcpy(&(array_out[row * max_column_count]),
|
|
&(array_out[(size_t)i * max_column_count]),
|
|
((size_t)*n + 1) * sizeof(LDBLE));
|
|
}
|
|
row_back_l[row] = i;
|
|
row++;
|
|
k1++;
|
|
}
|
|
|
|
/*
|
|
* Eliminate unnecessary equality eqns
|
|
*/
|
|
l1 = 0;
|
|
for (i = *k; i < (*k + *l); i++)
|
|
{
|
|
nonzero = FALSE;
|
|
for (j = 0; j < *n; j++)
|
|
{
|
|
if (equal(array_out[i * max_column_count + j], 0.0, toler) ==
|
|
FALSE)
|
|
{
|
|
nonzero = TRUE;
|
|
break;
|
|
}
|
|
}
|
|
if (nonzero == FALSE)
|
|
continue;
|
|
/*
|
|
if (memcmp(&(array_out[i * max_column_count]), &(zero[0]),
|
|
(size_t) (*n) * sizeof(LDBLE)) == 0) {
|
|
continue;
|
|
}
|
|
*/
|
|
/*
|
|
memcpy(&(array_out[row * max_column_count]), &(array_out[i * max_column_count]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
*/
|
|
|
|
if (i > row)
|
|
{
|
|
if ((row*max_column_count + (*n + 1)) >= (i * max_column_count))
|
|
{
|
|
assert(false);
|
|
}
|
|
memcpy(&(array_out[row * max_column_count]),
|
|
&(array_out[(size_t)i * max_column_count]),
|
|
((size_t)*n + 1) * sizeof(LDBLE));
|
|
}
|
|
row_back_l[row] = i;
|
|
row++;
|
|
l1++;
|
|
}
|
|
/*
|
|
* Eliminate unnecessary inequality eqns
|
|
*/
|
|
m1 = 0;
|
|
for (i = (*k + *l); i < (*k + *l + *m); i++)
|
|
{
|
|
nonzero = FALSE;
|
|
for (j = 0; j < *n; j++)
|
|
{
|
|
if (equal(array_out[i * max_column_count + j], 0.0, toler) ==
|
|
FALSE)
|
|
{
|
|
nonzero = TRUE;
|
|
break;
|
|
}
|
|
}
|
|
if (nonzero == FALSE)
|
|
continue;
|
|
/*
|
|
if (memcmp(&(array_out[i * max_column_count]), &(zero[0]),
|
|
(size_t) (*n) * sizeof(LDBLE)) == 0) {
|
|
continue;
|
|
}
|
|
*/
|
|
/*
|
|
memcpy(&(array_out[row * max_column_count]), &(array_out[i * max_column_count]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
*/
|
|
|
|
if (i > row)
|
|
{
|
|
if ((row*max_column_count + (*n + 1)) >= (i * max_column_count))
|
|
{
|
|
assert(false);
|
|
}
|
|
memcpy(&(array_out[(size_t)row * max_column_count]),
|
|
&(array_out[(size_t)i * max_column_count]),
|
|
((size_t)*n + 1) * sizeof(LDBLE));
|
|
}
|
|
row_back_l[row] = i;
|
|
row++;
|
|
m1++;
|
|
}
|
|
|
|
*k = k1;
|
|
*l = l1;
|
|
*m = m1;
|
|
/*
|
|
* Scale all inequality rows
|
|
*/
|
|
|
|
for (i = *k + *l; i < *k + *l + *m; i++)
|
|
{
|
|
for (j = 0; j < *n + 1; j++)
|
|
{
|
|
array_out[i * max_column_count + j] *= SCALE_ALL;
|
|
}
|
|
}
|
|
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
check_solns(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Check_solns checks that each solution can be charge balanced within
|
|
* the given constraints. If not, it is an error and the program will
|
|
* terminate.
|
|
*/
|
|
int i;
|
|
size_t j;
|
|
int k, l, m, n;
|
|
int return_value;
|
|
unsigned long bits;
|
|
LDBLE error2;
|
|
|
|
memcpy((void *) &(min_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(max_delta[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
|
|
/*
|
|
* Switch bits so that phases are high and solutions are low
|
|
*/
|
|
return_value = OK;
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
bits = 0;
|
|
bits += 1 << (inv_ptr->phases.size() + i);
|
|
/*
|
|
* Check for feasibility of charge balance with given uncertainties
|
|
*/
|
|
k = (int)row_mb; /* rows in A */
|
|
l = (int)(row_epsilon - row_mb); /* rows in C */
|
|
m = (int)(count_rows - row_epsilon); /* rows in E */
|
|
n = (int)count_unknowns; /* number of variables */
|
|
/* debug
|
|
output_msg(sformatf( "\nColumns\n"));
|
|
for (j = 0; j < n; j++) {
|
|
output_msg(sformatf( "\t%d\t%s\n", j, col_name[j]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nRows\n"));
|
|
for (j = 0; j < k + l + m; j++) {
|
|
output_msg(sformatf( "\t%d\t%s\n", j, row_name[j]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nA and B arrays:\n\n"));
|
|
array_print(array, k + l + m,
|
|
n + 1, max_column_count);
|
|
*/
|
|
/*
|
|
* Copy equations
|
|
*/
|
|
memcpy((void *) &(array1[0]), (void *) &(my_array[0]),
|
|
(size_t) max_column_count * max_row_count * sizeof(LDBLE));
|
|
memcpy((void *) &(delta2[0]), (void *) &(delta[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
memcpy((void *) &(inv_res[0]), (void *) &(inv_zero[0]),
|
|
(size_t) max_row_count * sizeof(LDBLE));
|
|
|
|
/*
|
|
* Keep optimization
|
|
*/
|
|
/*
|
|
* Zero out mass balance rows and fraction rows
|
|
*/
|
|
for (j = row_mb; j < row_charge; j++)
|
|
{
|
|
memcpy((void *) &(array1[j * max_column_count]),
|
|
(void *) &(inv_zero[0]),
|
|
max_column_count * sizeof(LDBLE));
|
|
}
|
|
/*
|
|
* Set fraction of solution to 1.0
|
|
*/
|
|
array1[(row_charge - 1) * max_column_count + i] = 1.0;
|
|
array1[(row_charge - 1) * max_column_count + n] = 1.0;
|
|
|
|
/*
|
|
* Zero out charge balance rows for other solutions
|
|
*/
|
|
for (j = 0; j < inv_ptr->count_solns; j++)
|
|
{
|
|
if (j == i)
|
|
continue;
|
|
memcpy((void *) &(array1[(row_charge + j) * max_column_count]),
|
|
(void *) &(inv_zero[0]),
|
|
(size_t) max_column_count * sizeof(LDBLE));
|
|
}
|
|
|
|
/*
|
|
* Zero out isotope mole balance
|
|
*/
|
|
for (size_t j = row_isotopes; j < row_epsilon; j++)
|
|
{
|
|
memcpy((void *) &(array1[j * max_column_count]),
|
|
(void *) &(inv_zero[0]),
|
|
max_column_count * sizeof(LDBLE));
|
|
}
|
|
|
|
/*
|
|
* Zero out isotope uncertainties
|
|
*/
|
|
for (size_t j = row_isotope_epsilon; j < count_rows; j++)
|
|
{
|
|
memcpy((void *) &(array1[j * max_column_count]),
|
|
(void *) &(inv_zero[0]),
|
|
max_column_count * sizeof(LDBLE));
|
|
}
|
|
/*
|
|
* Can`t Zero out epsilon constraint rows for other solutions because not sure which
|
|
* are which
|
|
*/
|
|
|
|
shrink(inv_ptr, &array1[0], &array1[0],
|
|
&k, &l, &m, &n, bits, &delta2[0], &col_back[0], &row_back[0]);
|
|
/* Debug
|
|
|
|
output_msg(sformatf( "\nColumns\n"));
|
|
for (j = 0; j < n; j++) {
|
|
output_msg(sformatf( "\t%d\t%s\n", j, col_name[col_back[j]]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nRows\n"));
|
|
for (j = 0; j < k + l + m; j++) {
|
|
output_msg(sformatf( "\t%d\t%s\n", j, row_name[row_back[j]]));
|
|
}
|
|
|
|
output_msg(sformatf( "\nA and B arrays:\n\n"));
|
|
array_print(array1, k + l + m,
|
|
n + 1, max_column_count);
|
|
|
|
output_msg(sformatf( "\nInput delta vector:\n"));
|
|
for (j=0; j < n; j++) {
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", j, col_name[col_back[j]], delta2[j]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
*/
|
|
|
|
kode = 1;
|
|
iter = 200;
|
|
count_calls++;
|
|
cl1(k, l, m, n, (int)nklmd, (int)n2d, &array1[0], &kode, toler, &iter,
|
|
&delta2[0], &inv_res[0], &error2, &inv_cu[0], &inv_iu[0], &inv_is[0], TRUE);
|
|
|
|
if (kode != 0)
|
|
{
|
|
error_string = sformatf(
|
|
"Not possible to balance solution %d with input uncertainties.",
|
|
inv_ptr->solns[i]);
|
|
error_msg(error_string, CONTINUE);
|
|
return_value = ERROR;
|
|
}
|
|
|
|
/* Debug
|
|
output_msg(sformatf( "kode: %d\titer: %d\terror: %e\n", kode, iter, error));
|
|
output_msg(sformatf( "k, l, m, n: %d\t%d\t%d\t%d\n", k, l, m, n));
|
|
|
|
output_msg(sformatf( "\nsolution vector %s\n", col_name[i]));
|
|
for (j = 0; j < n; j++) {
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e", j, col_name[col_back[j]], delta2[j]));
|
|
output_msg(sformatf( "\n"));
|
|
}
|
|
|
|
output_msg(sformatf( "\nresidual vector:\n"));
|
|
for (j = 0; j < (k + l + m); j++) {
|
|
output_msg(sformatf( "%6d %-12.12s %10.2e\n", j, row_name[row_back[j]], inv_res[j]));
|
|
}
|
|
*/
|
|
}
|
|
return (return_value);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
post_mortem(void)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Post_mortem simply identifies which equality and inequality of the
|
|
* array have not been satisfied.
|
|
*
|
|
*/
|
|
LDBLE sum;
|
|
/*
|
|
* Check equalities
|
|
*/
|
|
output_msg(sformatf(
|
|
"\nPost_mortem examination of inverse modeling:\n\n"));
|
|
for (size_t i = row_mb; i < row_epsilon; i++)
|
|
{
|
|
sum = 0;
|
|
for (size_t j = 0; j < count_unknowns; j++)
|
|
{
|
|
sum += inv_delta1[j] * my_array[i * max_column_count + j];
|
|
}
|
|
|
|
if (equal(sum, my_array[(i * max_column_count) + count_unknowns], toler)
|
|
== FALSE)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: equality not satisfied for %s, %e.\n",
|
|
row_name[i],
|
|
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
|
|
}
|
|
}
|
|
/*
|
|
* Check inequalities
|
|
*/
|
|
for (size_t i = row_epsilon; i < count_rows; i++)
|
|
{
|
|
sum = 0;
|
|
for (size_t j = 0; j < count_unknowns; j++)
|
|
{
|
|
sum += inv_delta1[j] * my_array[i * max_column_count + j];
|
|
}
|
|
|
|
if (sum > my_array[(i * max_column_count) + count_unknowns] + toler)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: inequality not satisfied for %s, %e\n",
|
|
row_name[i],
|
|
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
|
|
}
|
|
}
|
|
/*
|
|
* Check dissolution/precipitation constraints
|
|
*/
|
|
for (size_t i = 0; i < count_unknowns; i++)
|
|
{
|
|
if (delta_save[i] > 0.5 && inv_delta1[i] < -toler)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: Dissolution/precipitation constraint not satisfied for column %d, %s, %e.\n",
|
|
i, col_name[i], (double) inv_delta1[i]));
|
|
}
|
|
else if (delta_save[i] < -0.5 && inv_delta1[i] > toler)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: Dissolution/precipitation constraint not satisfied for column %d, %s, %e.\n",
|
|
i, col_name[i], (double) inv_delta1[i]));
|
|
}
|
|
}
|
|
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
bool Phreeqc::
|
|
test_cl1_solution(void)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* checks that equality and inequalities are satisfied
|
|
*
|
|
*/
|
|
int i;
|
|
LDBLE sum;
|
|
/*
|
|
* Check equalities
|
|
*/
|
|
bool rv = true;
|
|
if (debug_inverse)
|
|
{
|
|
output_msg(sformatf(
|
|
"\nTesting cl1 inverse modeling:\n\n"));
|
|
}
|
|
for (size_t i = row_mb; i < row_epsilon; i++)
|
|
{
|
|
sum = 0;
|
|
for (size_t j = 0; j < count_unknowns; j++)
|
|
{
|
|
sum += inv_delta1[j] * my_array[i * max_column_count + j];
|
|
}
|
|
|
|
if (equal(sum, my_array[(i * max_column_count) + count_unknowns], toler) == FALSE)
|
|
{
|
|
if (debug_inverse)
|
|
{
|
|
output_msg(sformatf("\tERROR: equality not satisfied for %s, %e.\n", row_name[i],
|
|
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
|
|
}
|
|
rv = false;
|
|
}
|
|
}
|
|
/*
|
|
* Check inequalities
|
|
*/
|
|
for (size_t i = row_epsilon; i < count_rows; i++)
|
|
{
|
|
sum = 0;
|
|
for (size_t j = 0; j < count_unknowns; j++)
|
|
{
|
|
sum += inv_delta1[j] * my_array[i * max_column_count + j];
|
|
}
|
|
|
|
if (sum > my_array[(i * max_column_count) + count_unknowns] + toler)
|
|
{
|
|
if (debug_inverse)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: inequality not satisfied for %s, %e\n",
|
|
row_name[i],
|
|
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
|
|
}
|
|
rv = false;
|
|
}
|
|
}
|
|
/*
|
|
* Check dissolution/precipitation constraints
|
|
*/
|
|
for (i = 0; i < count_unknowns; i++)
|
|
{
|
|
if (delta_save[i] > 0.5 && inv_delta1[i] < -toler)
|
|
{
|
|
if (debug_inverse)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: Dissolution/precipitation constraint not satisfied for column %d, %s, %e.\n",
|
|
i, col_name[i], (double) inv_delta1[i]));
|
|
}
|
|
rv = false;
|
|
|
|
}
|
|
else if (delta_save[i] < -0.5 && inv_delta1[i] > toler)
|
|
{
|
|
if (debug_inverse)
|
|
{
|
|
output_msg(sformatf(
|
|
"\tERROR: Dissolution/precipitation constraint not satisfied for column %d, %s, %e.\n",
|
|
i, col_name[i], (double) inv_delta1[i]));
|
|
}
|
|
rv = false;
|
|
}
|
|
}
|
|
|
|
return (rv);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
carbon_derivs(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
int i, j, temp;
|
|
LDBLE c_uncertainty, d_carbon, alk_plus, alk_minus;
|
|
cxxSolution *solution_ptr_orig, *solution_ptr;
|
|
|
|
inv_ptr->dalk_dph.resize(inv_ptr->count_solns);
|
|
inv_ptr->dalk_dc.resize(inv_ptr->count_solns);
|
|
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
solution_ptr_orig = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
if (solution_ptr_orig == NULL)
|
|
{
|
|
error_string = sformatf( "Solution %d for inverse "
|
|
"modeling not found.", inv_ptr->solns[i]);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
/*
|
|
* Find carbon uncertainty
|
|
*/
|
|
c_uncertainty = 0;
|
|
d_carbon = 0;
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
if (inv_ptr->elts[j].master == s_co3->secondary)
|
|
{
|
|
c_uncertainty = inv_ptr->elts[j].uncertainties[i];
|
|
break;
|
|
}
|
|
}
|
|
if (c_uncertainty < 0.0)
|
|
{
|
|
d_carbon = -c_uncertainty;
|
|
}
|
|
else if (c_uncertainty > 0.0)
|
|
{
|
|
cxxNameDouble::iterator kit = solution_ptr_orig->Get_totals().begin();
|
|
for ( ; kit != solution_ptr_orig->Get_totals().end(); kit++)
|
|
{
|
|
if (strcmp(kit->first.c_str(), "C(4)") == 0)
|
|
{
|
|
d_carbon = kit->second /
|
|
solution_ptr_orig->Get_mass_water() * c_uncertainty;
|
|
break;
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Make four copies of solution
|
|
* Modify ph and carbon in solutions
|
|
*/
|
|
set_ph_c(inv_ptr, i, solution_ptr_orig, -5, 0.0, 1.0, 0.0);
|
|
set_ph_c(inv_ptr, i, solution_ptr_orig, -4, 0.0, -1.0, 0.0);
|
|
if (c_uncertainty != 0)
|
|
{
|
|
set_ph_c(inv_ptr, i, solution_ptr_orig, -3, d_carbon, 0.0, 1.0);
|
|
set_ph_c(inv_ptr, i, solution_ptr_orig, -2, d_carbon, 0.0, -1.0);
|
|
}
|
|
/* */
|
|
temp = pr.all;
|
|
pr.all = FALSE;
|
|
initial_solutions(FALSE);
|
|
pr.all = temp;
|
|
/*
|
|
* dAlk/dpH
|
|
*/
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -5);
|
|
alk_plus = solution_ptr->Get_total_alkalinity();
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -4);
|
|
alk_minus = solution_ptr->Get_total_alkalinity();
|
|
inv_ptr->dalk_dph[i] = (alk_plus - alk_minus) /
|
|
(2.0 * inv_ptr->ph_uncertainties[i]);
|
|
/*
|
|
* dAlk/dC
|
|
*/
|
|
if (d_carbon != 0)
|
|
{
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -3);
|
|
alk_plus = solution_ptr->Get_total_alkalinity();
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -2);
|
|
alk_minus = solution_ptr->Get_total_alkalinity();
|
|
inv_ptr->dalk_dc[i] = (alk_plus - alk_minus) / (2.0 * d_carbon);
|
|
}
|
|
else
|
|
{
|
|
inv_ptr->dalk_dc[i] = 0.0;
|
|
}
|
|
if (debug_inverse == TRUE)
|
|
{
|
|
output_msg(sformatf( "dAlk/dph = %e\tdAlk/dC = %e\n",
|
|
(double) inv_ptr->dalk_dph[i],
|
|
(double) inv_ptr->dalk_dc[i]));
|
|
}
|
|
}
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
set_ph_c(class inverse *inv_ptr,
|
|
int i,
|
|
cxxSolution *solution_ptr_orig,
|
|
int n_user_new, LDBLE d_carbon, LDBLE ph_factor, LDBLE c_factor)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
int n_user_orig;
|
|
cxxSolution *solution_ptr;
|
|
|
|
n_user_orig = inv_ptr->solns[i];
|
|
Utilities::Rxn_copy(Rxn_solution_map, n_user_orig, n_user_new);
|
|
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, n_user_new);
|
|
solution_ptr->Set_new_def(true);
|
|
solution_ptr->Create_initial_data();
|
|
solution_ptr->Set_n_user_end(n_user_new);
|
|
LDBLE ph = solution_ptr->Get_ph();
|
|
ph += inv_ptr->ph_uncertainties[i] * ph_factor;
|
|
solution_ptr->Set_ph(ph);
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
cxxISolutionComp temp_comp;
|
|
temp_comp.Set_description(jit->first.c_str());
|
|
temp_comp.Set_input_conc(jit->second / solution_ptr_orig->Get_mass_water());
|
|
temp_comp.Set_units("Mol/kgw");
|
|
if (strcmp(jit->first.c_str(), "C(4)") == 0)
|
|
{
|
|
LDBLE c = temp_comp.Get_input_conc();
|
|
c += d_carbon * c_factor;
|
|
temp_comp.Set_input_conc(c);
|
|
}
|
|
solution_ptr->Get_initial_data()->Get_comps()[jit->first] = temp_comp;
|
|
}
|
|
solution_ptr->Get_totals().clear();
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
isotope_balance_equation(class inverse *inv_ptr, int row, int n)
|
|
/* ---------------------------------------------------------------------- */
|
|
/*
|
|
* routine fills in an isotope balance equation
|
|
*
|
|
* row is the row in array that needs to be filled
|
|
* n is the isotope number in inv_ptr
|
|
*/
|
|
{
|
|
int i, j, k;
|
|
LDBLE isotope_number;
|
|
size_t column;
|
|
LDBLE f;
|
|
class master *primary_ptr;
|
|
cxxSolution *solution_ptr;
|
|
/*
|
|
* Determine primary master species and isotope number for
|
|
* isotope mass-balance equation
|
|
*/
|
|
column = 0;
|
|
primary_ptr = master_bsearch_primary(inv_ptr->isotopes[n].elt_name);
|
|
isotope_number = inv_ptr->isotopes[n].isotope_number;
|
|
/* isotope element must be defined */
|
|
if (primary_ptr == NULL)
|
|
{
|
|
error_string = sformatf(
|
|
"In isotope calculation: element not defined: %s.",
|
|
inv_ptr->isotopes[n].elt_name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
}
|
|
|
|
/* isotope element must be primary */
|
|
if (primary_ptr->primary != TRUE)
|
|
{
|
|
error_string = sformatf( "Isotope mass-balance may only be used"
|
|
" for total element concentrations.\n"
|
|
"Secondary species not allowed: %s.",
|
|
inv_ptr->isotopes[n].elt_name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
}
|
|
|
|
/*
|
|
* Fill in terms for each solution
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (i == (inv_ptr->count_solns - 1))
|
|
{
|
|
f = -1.0;
|
|
}
|
|
else
|
|
{
|
|
f = 1.0;
|
|
}
|
|
|
|
/* mixing fraction term */
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
std::map < std::string, cxxSolutionIsotope >::iterator jit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; jit != solution_ptr->Get_isotopes().end(); jit++)
|
|
{
|
|
class master *primary_jit = master_bsearch_primary(jit->second.Get_elt_name().c_str());
|
|
if (primary_jit == primary_ptr &&
|
|
jit->second.Get_isotope_number() == isotope_number)
|
|
{
|
|
my_array[(size_t)row * max_column_count + (size_t)i] +=
|
|
f * jit->second.Get_total() * jit->second.Get_ratio();
|
|
}
|
|
}
|
|
|
|
/* epsilon of total moles of element valence * ratio */
|
|
jit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; jit != solution_ptr->Get_isotopes().end(); jit++)
|
|
{
|
|
/* What to do with H and O, skip for now ??? */
|
|
if (primary_ptr == s_hplus->primary
|
|
|| primary_ptr == s_h2o->primary)
|
|
continue;
|
|
class master *master_jit = master_bsearch(jit->second.Get_elt_name().c_str());
|
|
class master *primary_jit = master_bsearch_primary(jit->second.Get_elt_name().c_str());
|
|
if (primary_jit == primary_ptr &&
|
|
jit->second.Get_isotope_number() == isotope_number)
|
|
{
|
|
/* find column of master for solution i */
|
|
for (k = 0; k < inv_ptr->elts.size(); k++)
|
|
{
|
|
if (master_jit == inv_ptr->elts[k].master)
|
|
break;
|
|
}
|
|
column = col_epsilon + (k * inv_ptr->count_solns) + i;
|
|
my_array[(size_t)row * max_column_count + (size_t)column] +=
|
|
f * jit->second.Get_ratio();
|
|
}
|
|
}
|
|
|
|
/* epsilon of ratio * total of element valence */
|
|
jit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; jit != solution_ptr->Get_isotopes().end(); jit++)
|
|
{
|
|
class master *master_jit = master_bsearch(jit->second.Get_elt_name().c_str());
|
|
class master *primary_jit = master_bsearch_primary(jit->second.Get_elt_name().c_str());
|
|
if (primary_jit == primary_ptr &&
|
|
jit->second.Get_isotope_number() == isotope_number)
|
|
{
|
|
|
|
/* find column of epsilon for ratio of valence */
|
|
for (k = 0; k < inv_ptr->isotope_unknowns.size(); k++)
|
|
{
|
|
if (master_jit ==
|
|
inv_ptr->isotope_unknowns[k].master
|
|
&& jit->second.Get_isotope_number() ==
|
|
inv_ptr->isotope_unknowns[k].isotope_number)
|
|
{
|
|
column =
|
|
col_isotopes +
|
|
(i * inv_ptr->isotope_unknowns.size()) + k;
|
|
}
|
|
}
|
|
my_array[(size_t)row * max_column_count + (size_t)column] +=
|
|
f * jit->second.Get_total();
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Fill in terms for each phase
|
|
*/
|
|
for (i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (inv_ptr->phases[i].isotopes.size() == 0)
|
|
continue;
|
|
std::vector<class isotope>& isotope_ref = inv_ptr->phases[i].isotopes;
|
|
for (j = 0; j < inv_ptr->phases[i].isotopes.size(); j++)
|
|
{
|
|
if (isotope_ref[j].primary == primary_ptr &&
|
|
isotope_ref[j].isotope_number == isotope_number)
|
|
{
|
|
/* term for alpha phase unknowns */
|
|
column = col_phases + i;
|
|
my_array[(size_t)row * max_column_count + (size_t)column] =
|
|
isotope_ref[j].ratio * isotope_ref[j].coef;
|
|
/* term for phase isotope uncertainty unknown */
|
|
column = col_phase_isotopes + i * inv_ptr->isotopes.size() + (size_t)n;
|
|
my_array[(size_t)row * max_column_count + column] = isotope_ref[j].coef;
|
|
break;
|
|
}
|
|
}
|
|
|
|
}
|
|
return OK;
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
bool Phreeqc::
|
|
set_isotope_unknowns(class inverse* inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Go through elements for which isotope balances are requested
|
|
* and make a array of isotope structures
|
|
* return total number of isotope unknowns and structure array
|
|
*/
|
|
int i, k;
|
|
LDBLE isotope_number;
|
|
class master* primary_ptr;
|
|
size_t count_isotopes;
|
|
std::vector<class isotope>& isotopes = inv_ptr->isotope_unknowns;
|
|
|
|
if (inv_ptr->isotopes.size() == 0)
|
|
{
|
|
isotopes.clear();
|
|
return true;
|
|
}
|
|
count_isotopes = 0;
|
|
|
|
for (i = 0; i < inv_ptr->isotopes.size(); i++)
|
|
{
|
|
primary_ptr = master_bsearch(inv_ptr->isotopes[i].elt_name);
|
|
isotope_number = inv_ptr->isotopes[i].isotope_number;
|
|
if (primary_ptr == NULL)
|
|
{
|
|
error_string = sformatf(
|
|
"Element not found for isotope calculation: %s.",
|
|
inv_ptr->isotopes[i].elt_name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
break;
|
|
}
|
|
if (primary_ptr->primary != TRUE)
|
|
{
|
|
error_string = sformatf("Isotope mass-balance may only be used"
|
|
" for total element concentrations.\n"
|
|
"Secondary species not allowed: %s.",
|
|
inv_ptr->isotopes[i].elt_name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
break;
|
|
}
|
|
|
|
/* nonredox element */
|
|
if (primary_ptr->s->secondary == NULL)
|
|
{
|
|
isotopes.resize(count_isotopes + 1);
|
|
isotopes[count_isotopes].primary = primary_ptr;
|
|
isotopes[count_isotopes].master = primary_ptr;
|
|
isotopes[count_isotopes].isotope_number = isotope_number;
|
|
isotopes[count_isotopes].elt_name = primary_ptr->elt->name;
|
|
count_isotopes++;
|
|
|
|
/* redox element */
|
|
}
|
|
else
|
|
{
|
|
|
|
/* find master */
|
|
for (k = 0; k < (int)master.size(); k++)
|
|
{
|
|
if (master[k] == primary_ptr)
|
|
break;
|
|
}
|
|
|
|
/* sum all secondary for master */
|
|
k++;
|
|
for (; k < (int)master.size(); k++)
|
|
{
|
|
isotopes.resize(count_isotopes + 1);
|
|
isotopes[count_isotopes].primary = primary_ptr;
|
|
isotopes[count_isotopes].master = master[k];
|
|
isotopes[count_isotopes].isotope_number = isotope_number;
|
|
isotopes[count_isotopes].elt_name = master[k]->elt->name;
|
|
count_isotopes++;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
check_isotopes(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Go through elements for which isotope balances are requested
|
|
* and make sure each solution has isotope ratios defined
|
|
*/
|
|
int i, ii, j, k, l;
|
|
int err, found_isotope;
|
|
LDBLE isotope_number;
|
|
class master *master_ptr, *primary_ptr;
|
|
cxxSolution *solution_ptr;
|
|
class phase *phase_ptr;
|
|
char token[MAX_LENGTH];
|
|
|
|
/*
|
|
* Check solutions for necessary isotope data
|
|
*/
|
|
for (j = 0; j < inv_ptr->count_solns; j++)
|
|
{
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[j]);
|
|
xsolution_zero();
|
|
add_solution(solution_ptr, 1.0, 1.0);
|
|
/*
|
|
* Go through inverse isotopes and make sure isotope data for each solution
|
|
* inv_ptr->isotopes has elements; inv_ptr->i_u has redox states and uncertainties
|
|
*/
|
|
for (i = 0; i < inv_ptr->isotopes.size(); i++)
|
|
{
|
|
err = FALSE;
|
|
primary_ptr = master_bsearch(inv_ptr->isotopes[i].elt_name);
|
|
isotope_number = inv_ptr->isotopes[i].isotope_number;
|
|
found_isotope = FALSE;
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
class master *primary_kit = master_bsearch_primary(kit->second.Get_elt_name().c_str());
|
|
if (primary_kit == primary_ptr &&
|
|
kit->second.Get_isotope_number() ==
|
|
isotope_number)
|
|
{
|
|
found_isotope = TRUE;
|
|
break;
|
|
}
|
|
}
|
|
if (found_isotope == TRUE)
|
|
continue;
|
|
|
|
/* did not find isotope, which is ok if element not in solution */
|
|
if (primary_ptr == s_h2o->primary
|
|
|| primary_ptr == s_hplus->primary)
|
|
{
|
|
err = TRUE;
|
|
}
|
|
else if (primary_ptr->total > 0)
|
|
{
|
|
err = TRUE;
|
|
}
|
|
if (err == TRUE)
|
|
{
|
|
error_string = sformatf(
|
|
"In solution %d, isotope ratio(s) are needed for element: %g%s.",
|
|
solution_ptr->Get_n_user(), (double) isotope_number,
|
|
primary_ptr->elt->name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
continue;
|
|
}
|
|
}
|
|
/*
|
|
* Go through solution isotopes and set uncertainties
|
|
*/
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
class master *master_kit = master_bsearch(kit->second.Get_elt_name().c_str());
|
|
class master *primary_kit = master_bsearch_primary(kit->second.Get_elt_name().c_str());
|
|
kit->second.Set_x_ratio_uncertainty(NAN);
|
|
/*
|
|
* Search for secondary or primary master in inverse uncertainties
|
|
*/
|
|
ii = -1;
|
|
for (i = 0; i < inv_ptr->i_u.size(); i++)
|
|
{
|
|
master_ptr = master_bsearch(inv_ptr->i_u[i].elt_name);
|
|
if (master_ptr == master_kit)
|
|
{
|
|
ii = i;
|
|
break;
|
|
}
|
|
if (master_ptr == primary_kit)
|
|
{
|
|
ii = i;
|
|
}
|
|
}
|
|
/* solution isotope data not being used in inverse */
|
|
if (ii == -1)
|
|
continue;
|
|
|
|
i = ii;
|
|
/* use inverse-defined uncertainties first */
|
|
#ifdef NPP
|
|
if (j < inv_ptr->i_u[i].uncertainties.size()
|
|
&& !isnan(inv_ptr->i_u[i].uncertainties[j]))
|
|
#else
|
|
if (j < inv_ptr->i_u[i].uncertainties.size()
|
|
&& inv_ptr->i_u[i].uncertainties[j] != NAN)
|
|
#endif
|
|
{
|
|
kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]);
|
|
|
|
/* use solution-defined uncertainties second */
|
|
}
|
|
#ifdef NPP
|
|
else if (inv_ptr->i_u[i].uncertainties.size() > 0
|
|
&& !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
|
|
#else
|
|
else if (inv_ptr->i_u[i].uncertainties.size() > 0
|
|
&& inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN)
|
|
#endif
|
|
{
|
|
kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]);
|
|
|
|
/* use solution-defined uncertainties second */
|
|
}
|
|
#ifdef NPP
|
|
else if (!isnan(kit->second.Get_ratio_uncertainty()))
|
|
#else
|
|
else if (kit->second.Get_ratio_uncertainty() != NAN)
|
|
#endif
|
|
{
|
|
kit->second.Set_x_ratio_uncertainty(
|
|
kit->second.Get_ratio_uncertainty());
|
|
/* use isotope defaults third */
|
|
}
|
|
else
|
|
{
|
|
sprintf(token, "%g%s",
|
|
(double) kit->second.Get_isotope_number(),
|
|
kit->second.Get_elt_name().c_str());
|
|
for (l = 0; l < count_iso_defaults; l++)
|
|
{
|
|
if (strcmp(token, iso_defaults[l].name) == 0)
|
|
{
|
|
kit->second.Set_x_ratio_uncertainty(
|
|
iso_defaults[l].uncertainty);
|
|
error_string = sformatf(
|
|
"Solution %d, element %g%s: default isotope ratio uncertainty is used, %g.",
|
|
solution_ptr->Get_n_user(),
|
|
(double) kit->second.Get_isotope_number(),
|
|
kit->second.Get_elt_name().c_str(),
|
|
kit->second.Get_x_ratio_uncertainty());
|
|
warning_msg(error_string);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
#ifdef NPP
|
|
if (isnan(kit->second.Get_x_ratio_uncertainty()))
|
|
#else
|
|
if (kit->second.Get_x_ratio_uncertainty() == NAN)
|
|
#endif
|
|
{
|
|
error_string = sformatf(
|
|
"In solution %d, isotope ratio uncertainty is needed for element: %g%s.",
|
|
solution_ptr->Get_n_user(),
|
|
(double) kit->second.Get_isotope_number(),
|
|
kit->second.Get_elt_name().c_str());
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Check phases for necessary isotope data
|
|
*/
|
|
for (j = 0; j < inv_ptr->phases.size(); j++)
|
|
{
|
|
for (i = 0; i < inv_ptr->isotopes.size(); i++)
|
|
{
|
|
primary_ptr = master_bsearch(inv_ptr->isotopes[i].elt_name);
|
|
isotope_number = inv_ptr->isotopes[i].isotope_number;
|
|
found_isotope = FALSE;
|
|
for (k = 0; k < inv_ptr->phases[j].isotopes.size(); k++)
|
|
{
|
|
if (inv_ptr->phases[j].isotopes[k].primary == primary_ptr &&
|
|
inv_ptr->phases[j].isotopes[k].isotope_number ==
|
|
isotope_number)
|
|
{
|
|
found_isotope = TRUE;
|
|
break;
|
|
}
|
|
}
|
|
if (found_isotope == TRUE)
|
|
continue;
|
|
|
|
/* did not find isotope, which is ok if element not in solution */
|
|
phase_ptr = inv_ptr->phases[j].phase;
|
|
k = 0;
|
|
while (phase_ptr->next_elt[k].elt != NULL)
|
|
{
|
|
if (phase_ptr->next_elt[k].elt->primary == primary_ptr)
|
|
{
|
|
if (s_hplus->primary == primary_ptr ||
|
|
s_h2o->primary == primary_ptr)
|
|
{
|
|
k++;
|
|
continue;
|
|
}
|
|
else
|
|
{
|
|
error_string = sformatf(
|
|
"In phase %s, isotope ratio(s) are needed for element: %g%s.",
|
|
phase_ptr->name, (double) isotope_number,
|
|
primary_ptr->elt->name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
break;
|
|
}
|
|
}
|
|
k++;
|
|
}
|
|
}
|
|
}
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
phase_isotope_inequalities(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
size_t column;
|
|
char token[MAX_LENGTH];
|
|
if (inv_ptr->isotopes.size() <= 0)
|
|
return OK;
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
if (inv_ptr->phases[i].isotopes.size() == 0)
|
|
continue;
|
|
|
|
for (size_t j = 0; j < inv_ptr->phases[i].isotopes.size(); j++)
|
|
{
|
|
/* find index number */
|
|
size_t k = 0;
|
|
for (k = 0; k < inv_ptr->isotopes.size(); k++)
|
|
{
|
|
if (inv_ptr->phases[i].isotopes[j].elt_name ==
|
|
inv_ptr->isotopes[k].elt_name
|
|
&& inv_ptr->phases[i].isotopes[j].isotope_number ==
|
|
inv_ptr->isotopes[k].isotope_number)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
if (k >= inv_ptr->isotopes.size())
|
|
break;
|
|
column = col_phase_isotopes + i * inv_ptr->isotopes.size() + k;
|
|
/*
|
|
* zero column if uncertainty is zero
|
|
*/
|
|
if (inv_ptr->phases[i].isotopes[j].ratio_uncertainty == 0)
|
|
{
|
|
for (k = 0; k < count_rows; k++)
|
|
{
|
|
my_array[(size_t)k * max_column_count + (size_t)column] = 0.0;
|
|
}
|
|
continue;
|
|
}
|
|
|
|
/*
|
|
* optimization
|
|
*/
|
|
my_array[((size_t)column - (size_t)col_epsilon) * max_column_count + (size_t)column] =
|
|
SCALE_EPSILON / inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
|
|
/*
|
|
* two inequalities to account for absolute value
|
|
*/
|
|
/* for phases constrained to precipitate */
|
|
if (inv_ptr->phases[i].constraint == PRECIPITATE)
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)col_phases + (size_t)i] =
|
|
inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0;
|
|
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
|
|
"iso pos");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
my_array[count_rows * max_column_count + (size_t)col_phases + (size_t)i] =
|
|
inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
|
|
"iso neg");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* for phases constrained to dissolve */
|
|
}
|
|
else if (inv_ptr->phases[i].constraint == DISSOLVE)
|
|
{
|
|
my_array[count_rows * max_column_count + (size_t)col_phases + (size_t)i] =
|
|
-inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
|
|
my_array[count_rows * max_column_count + (size_t)column] = -1.0;
|
|
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
|
|
"iso pos");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
my_array[count_rows * max_column_count + (size_t)col_phases + (size_t)i] =
|
|
-inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
|
|
my_array[count_rows * max_column_count + (size_t)column] = 1.0;
|
|
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
|
|
"iso neg");
|
|
row_name[count_rows] = string_hsave(token);
|
|
count_rows++;
|
|
|
|
/* Error if phase is not constrained */
|
|
}
|
|
else
|
|
{
|
|
error_string = sformatf(
|
|
"In isotope calculations, all phases containing isotopes must be"
|
|
" constrained.\nPhase %s is not constrained.\n",
|
|
inv_ptr->phases[i].phase->name);
|
|
error_msg(error_string, CONTINUE);
|
|
input_error++;
|
|
continue;
|
|
}
|
|
}
|
|
}
|
|
return (OK);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
write_optimize_names(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
int i, j, row;
|
|
char token[MAX_LENGTH];
|
|
row = 0;
|
|
/*
|
|
* epsilons for analytical data
|
|
*/
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
sprintf(token, "%s %s %d", "optimize",
|
|
inv_ptr->elts[j].master->elt->name, inv_ptr->solns[i]);
|
|
row_name[row] = string_hsave(token);
|
|
row++;
|
|
}
|
|
}
|
|
/*
|
|
* pH
|
|
*/
|
|
if (carbon > 0)
|
|
{
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
sprintf(token, "%s %s %d", "optimize", "pH", inv_ptr->solns[i]);
|
|
row_name[row] = string_hsave(token);
|
|
row++;
|
|
}
|
|
}
|
|
/*
|
|
* water
|
|
*/
|
|
sprintf(token, "%s %s", "optimize", "water");
|
|
row_name[row] = string_hsave(token);
|
|
row++;
|
|
/*
|
|
* solution isotopes
|
|
*/
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
for (j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
sprintf(token, "%s %d%s %d", "optimize",
|
|
(int) inv_ptr->isotope_unknowns[j].isotope_number,
|
|
inv_ptr->isotope_unknowns[j].elt_name, inv_ptr->solns[i]);
|
|
row_name[row] = string_hsave(token);
|
|
row++;
|
|
}
|
|
}
|
|
/*
|
|
* phase isotopes
|
|
*/
|
|
|
|
for (i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
for (j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
sprintf(token, "%s %s %d%s", "optimize",
|
|
inv_ptr->phases[i].phase->name,
|
|
(int) inv_ptr->isotopes[j].isotope_number,
|
|
inv_ptr->isotopes[j].elt_name);
|
|
row_name[row] = string_hsave(token);
|
|
row++;
|
|
}
|
|
}
|
|
return OK;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
void Phreeqc::
|
|
dump_netpath(class inverse *inverse_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
std::string string;
|
|
//const char* cptr;
|
|
|
|
if (inverse_ptr->netpath == NULL)
|
|
return;
|
|
|
|
/* open file */
|
|
string = inverse_ptr->netpath;
|
|
if (replace(".lon", ".lon", string) != true)
|
|
{
|
|
string.append(".lon");
|
|
}
|
|
netpath_file = fopen(string.c_str(), "w");
|
|
if (netpath_file == NULL)
|
|
{
|
|
error_string = sformatf( "Can`t open file, %s.", inverse_ptr->netpath);
|
|
error_msg(error_string, STOP);
|
|
#if !defined(R_SO)
|
|
exit(4);
|
|
#endif
|
|
}
|
|
add_to_file("netpath.fil", inverse_ptr->netpath);
|
|
|
|
/* Header */
|
|
fprintf(netpath_file,
|
|
"2.14 # File format\n");
|
|
|
|
/* write out each solution */
|
|
std::map<int, cxxSolution>::iterator it = Rxn_solution_map.begin();
|
|
for ( ; it != Rxn_solution_map.end(); it++)
|
|
{
|
|
if (it->second.Get_n_user() < 0)
|
|
continue;
|
|
if (it->second.Get_description().size() > 0)
|
|
{
|
|
string = it->second.Get_description();
|
|
}
|
|
else
|
|
{
|
|
string = sformatf("Solution %d", it->second.Get_n_user());
|
|
}
|
|
fprintf(netpath_file, "4020%s\n", string.c_str());
|
|
//description = (char *) free_check_null(description);
|
|
/* lat/lon */
|
|
fprintf(netpath_file,
|
|
" # Lat/lon\n");
|
|
|
|
/* well number */
|
|
fprintf(netpath_file,
|
|
"%15d # Well number\n",
|
|
it->second.Get_n_user());
|
|
|
|
/* total number of wells */
|
|
fprintf(netpath_file,
|
|
"%15d # Total wells\n",
|
|
(int) Rxn_solution_map.size());
|
|
|
|
/* address */
|
|
fprintf(netpath_file,
|
|
" # Address1\n");
|
|
fprintf(netpath_file,
|
|
" # Address2\n");
|
|
fprintf(netpath_file,
|
|
" # Address3\n");
|
|
fprintf(netpath_file,
|
|
" # Address4\n");
|
|
fprintf(netpath_file,
|
|
" # Address5\n");
|
|
|
|
/* temperature */
|
|
fprintf(netpath_file,
|
|
"%15g # Temperature\n",
|
|
(double) it->second.Get_tc());
|
|
|
|
/* pH */
|
|
fprintf(netpath_file,
|
|
"%15g # pH\n",
|
|
(double) it->second.Get_ph());
|
|
|
|
/* DO */
|
|
print_total(netpath_file, &(it->second), "O(0)", "Dissolved Oxygen");
|
|
|
|
/* TDIC */
|
|
print_total(netpath_file, &(it->second), "C(4)", "TDIC");
|
|
|
|
/* Tritium */
|
|
print_isotope(netpath_file, &(it->second), "3H(1)", "Tritium");
|
|
|
|
/* H2S */
|
|
print_total(netpath_file, &(it->second), "S(-2)", "H2S");
|
|
|
|
/* Calcium */
|
|
print_total(netpath_file, &(it->second), "Ca", "Calcium");
|
|
|
|
/* Eh */
|
|
fprintf(netpath_file,
|
|
"%15g # Eh\n",
|
|
(double) (0.059 * it->second.Get_pe()));
|
|
|
|
/* Magnesium */
|
|
print_total(netpath_file, &(it->second), "Mg", "Magnesium");
|
|
|
|
/* Sodium */
|
|
print_total(netpath_file, &(it->second), "Na", "Sodium");
|
|
|
|
/* Potassium */
|
|
print_total(netpath_file, &(it->second), "K", "Potassium");
|
|
|
|
/* Chloride */
|
|
print_total(netpath_file, &(it->second), "Cl", "Chloride");
|
|
|
|
/* Sulfate */
|
|
print_total(netpath_file, &(it->second), "S(6)", "Sulfate");
|
|
|
|
/* Fluoride */
|
|
print_total(netpath_file, &(it->second), "F", "Fluoride");
|
|
|
|
/* Silica */
|
|
print_total(netpath_file, &(it->second), "Si", "Silica");
|
|
|
|
/* Bromide */
|
|
print_total(netpath_file, &(it->second), "Br", "Bromide");
|
|
|
|
/* Boron */
|
|
print_total(netpath_file, &(it->second), "B", "Boron");
|
|
|
|
/* Barium */
|
|
print_total(netpath_file, &(it->second), "Ba", "Barium");
|
|
|
|
/* Lithium */
|
|
print_total(netpath_file, &(it->second), "Li", "Lithium");
|
|
|
|
/* Strontium */
|
|
print_total(netpath_file, &(it->second), "Sr", "Strontium");
|
|
|
|
/* Iron */
|
|
print_total_multi(netpath_file, &(it->second), "Iron", "Fe", "Fe(2)",
|
|
"Fe(3)", "", "");
|
|
|
|
|
|
/* Manganese */
|
|
print_total_multi(netpath_file, &(it->second), "Manganese", "Mn",
|
|
"Mn(2)", "Mn(3)", "Mn(6)", "Mn(7)");
|
|
|
|
/* Nitrate */
|
|
print_total(netpath_file, &(it->second), "N(5)", "Nitrate");
|
|
|
|
/* Ammonium */
|
|
print_total_multi(netpath_file, &(it->second), "Ammonium", "N(-3)",
|
|
"Amm", "", "", "");
|
|
|
|
/* Phosphate */
|
|
print_total(netpath_file, &(it->second), "P", "Phosphate");
|
|
|
|
/* DOC */
|
|
print_total_multi(netpath_file, &(it->second), "DOC", "Fulvate",
|
|
"Humate", "", "", "");
|
|
|
|
/* Sp. Cond. */
|
|
fprintf(netpath_file,
|
|
" # Sp. Cond.\n");
|
|
|
|
/* Density */
|
|
fprintf(netpath_file,
|
|
" # Density\n");
|
|
|
|
/* Delta C-13 TDIC */
|
|
print_isotope(netpath_file, &(it->second), "13C(4)", "Delta C-13 TDIC");
|
|
|
|
/* C-14 TDIC */
|
|
print_isotope(netpath_file, &(it->second), "14C(4)", "C-14 TDIC");
|
|
|
|
/* Delta S-34 (SO4) */
|
|
print_isotope(netpath_file, &(it->second), "34S(6)",
|
|
"Delta S-34 (SO4)");
|
|
|
|
/* Delta S-34 (H2S) */
|
|
print_isotope(netpath_file, &(it->second), "34S(-2)",
|
|
"Delta S-34 (H2S)");
|
|
|
|
/* Delta Deuterium */
|
|
print_isotope(netpath_file, &(it->second), "2H(1)", "Delta Deuterium");
|
|
|
|
/* Delta O-18 */
|
|
print_isotope(netpath_file, &(it->second), "18O(-2)", "Delta O-18");
|
|
|
|
/* CH4 (aq) */
|
|
print_total(netpath_file, &(it->second), "C(-4)", "CH4 (aq)");
|
|
|
|
/* Sr 87/86 */
|
|
print_isotope(netpath_file, &(it->second), "87Sr", "Sr 87/86");
|
|
|
|
/* Al */
|
|
print_total(netpath_file, &(it->second), "Al", "Alumninum");
|
|
|
|
/* N2 (aq) */
|
|
print_total(netpath_file, &(it->second), "N(0)", "N2 (aq)");
|
|
|
|
/* N-15 of N2 (aq) */
|
|
print_isotope(netpath_file, &(it->second), "15N(0)", "N-15 of N2 (aq)");
|
|
|
|
/* N-15 of Nitrate */
|
|
print_isotope(netpath_file, &(it->second), "15N(5)", "N-15 of Nitrate");
|
|
|
|
/* N-15 of Ammonium */
|
|
print_isotope(netpath_file, &(it->second), "15N(-3)",
|
|
"N-15 of Ammonium");
|
|
|
|
/* Formation */
|
|
fprintf(netpath_file,
|
|
" # Formation\n");
|
|
|
|
}
|
|
if (netpath_file != NULL)
|
|
{
|
|
fclose(netpath_file);
|
|
netpath_file = NULL;
|
|
}
|
|
return;
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
LDBLE Phreeqc::
|
|
get_inv_total(cxxSolution *solution_ptr, const char *elt)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
if (strcmp(elt, jit->first.c_str()) == 0)
|
|
return jit->second;
|
|
}
|
|
return (0);
|
|
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
cxxSolutionIsotope *Phreeqc::get_isotope(cxxSolution *solution_ptr, const char *elt)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
std::string str_elt = elt;
|
|
std::map<std::string, cxxSolutionIsotope>::iterator it;
|
|
it = solution_ptr->Get_isotopes().find(str_elt);
|
|
if (it != solution_ptr->Get_isotopes().end())
|
|
{
|
|
return &(it->second);
|
|
}
|
|
return (NULL);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
void Phreeqc::
|
|
print_total(FILE * l_netpath_file, cxxSolution *solution_ptr,
|
|
const char *elt, const char *string)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
LDBLE moles = get_inv_total(solution_ptr, elt);
|
|
if (moles == 0)
|
|
{
|
|
fprintf(l_netpath_file,
|
|
" # %s\n",
|
|
string);
|
|
}
|
|
else
|
|
{
|
|
fprintf(l_netpath_file,
|
|
"%15g # %s\n",
|
|
(double) (1000 * moles / solution_ptr->Get_mass_water()),
|
|
string);
|
|
}
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
void Phreeqc::
|
|
print_isotope(FILE * l_netpath_file, cxxSolution *solution_ptr,
|
|
const char *elt, const char *string)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
cxxSolutionIsotope *iso_ptr;
|
|
iso_ptr = get_isotope(solution_ptr, elt);
|
|
if (iso_ptr == NULL)
|
|
{
|
|
fprintf(l_netpath_file,
|
|
" # %s\n",
|
|
string);
|
|
}
|
|
else
|
|
{
|
|
fprintf(l_netpath_file,
|
|
"%15g # %s\n",
|
|
(double) iso_ptr->Get_ratio(), string);
|
|
}
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
void Phreeqc::
|
|
print_total_multi(FILE * l_netpath_file, cxxSolution *solution_ptr,
|
|
const char *string, const char *elt0, const char *elt1,
|
|
const char *elt2, const char *elt3, const char *elt4)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
char elts[5][MAX_LENGTH];
|
|
LDBLE moles;
|
|
LDBLE sum;
|
|
int i, found;
|
|
|
|
strcpy(elts[0], elt0);
|
|
strcpy(elts[1], elt1);
|
|
strcpy(elts[2], elt2);
|
|
strcpy(elts[3], elt3);
|
|
strcpy(elts[4], elt4);
|
|
|
|
|
|
sum = 0;
|
|
found = FALSE;
|
|
for (i = 0; i < 5; i++)
|
|
{
|
|
moles = get_inv_total(solution_ptr, elts[i]);
|
|
if (moles == 0)
|
|
{
|
|
continue;
|
|
}
|
|
else
|
|
{
|
|
sum += moles;
|
|
found = TRUE;
|
|
}
|
|
}
|
|
if (found != TRUE)
|
|
{
|
|
fprintf(l_netpath_file,
|
|
" # %s\n",
|
|
string);
|
|
}
|
|
else
|
|
{
|
|
fprintf(l_netpath_file,
|
|
"%15g # %s\n",
|
|
(double) (1000 * sum / solution_ptr->Get_mass_water()), string);
|
|
}
|
|
return;
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
dump_netpath_pat(class inverse *inv_ptr)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
/*
|
|
* Prints model
|
|
*/
|
|
int i, j, k;
|
|
cxxSolution *solution_ptr, *solution_ptr_orig;
|
|
class master *master_ptr;
|
|
LDBLE d1, d2, d3;
|
|
LDBLE sum, sum1, sum_iso, d;
|
|
std::vector<double> array_save, l_delta_save;
|
|
size_t count_unknowns_save, max_row_count_save,
|
|
max_column_count_save, count_current_solutions;
|
|
int temp, temp_punch;
|
|
int solnmap[10][2];
|
|
FILE *model_file;
|
|
const class elt_list *next_elt;
|
|
int exch;
|
|
size_t column;
|
|
LDBLE f;
|
|
class rxn_token *rxn_ptr;
|
|
/*
|
|
* print solution data, epsilons, and revised data
|
|
*/
|
|
if (inv_ptr->pat == NULL)
|
|
return (OK);
|
|
|
|
array_save = my_array;
|
|
l_delta_save = delta;
|
|
count_unknowns_save = count_unknowns;
|
|
max_row_count_save = max_row_count;
|
|
max_column_count_save = max_column_count;
|
|
|
|
count_unknowns = 0;
|
|
max_row_count = 0;
|
|
max_column_count = 0;
|
|
|
|
count_current_solutions = 0;
|
|
count_inverse_models++;
|
|
|
|
for (i = 0; i < inv_ptr->count_solns; i++)
|
|
{
|
|
if (equal(inv_delta1[i], 0.0, toler) == TRUE)
|
|
continue;
|
|
solution_ptr_orig = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
|
|
Utilities::Rxn_copy(Rxn_solution_map, solution_ptr_orig->Get_n_user(), -6);
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -6);
|
|
xsolution_zero();
|
|
|
|
/* Adjust pH */
|
|
if (inv_ptr->carbon == TRUE)
|
|
{
|
|
d1 = solution_ptr->Get_ph();
|
|
d2 = inv_delta1[col_ph + i] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
solution_ptr->Set_ph(d3);
|
|
}
|
|
|
|
/* put original totals in master */
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
master_ptr = master_bsearch(jit->first.c_str());
|
|
master_ptr->total = jit->second;
|
|
}
|
|
|
|
/* ignore alkalinity */
|
|
/*master_alk->total = solution_ptr->total_alkalinity; */
|
|
|
|
/* update total in master */
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
if (inv_ptr->elts[j].master->s == s_eminus)
|
|
continue;
|
|
d1 = inv_ptr->elts[j].master->total;
|
|
d2 = inv_delta1[col_epsilon + j * inv_ptr->count_solns +
|
|
i] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
inv_ptr->elts[j].master->total = d3;
|
|
}
|
|
|
|
/* put updated total back in solution */
|
|
cxxNameDouble nd;
|
|
jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
master_ptr = master_bsearch(jit->first.c_str());
|
|
nd[jit->first] = master_ptr->total;
|
|
}
|
|
solution_ptr->Set_totals(nd);
|
|
|
|
/* update isotopes in solution */
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
{
|
|
/* adjustments to solution isotope composition */
|
|
for (j = 0; j < inv_ptr->isotope_unknowns.size(); j++)
|
|
{
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (inv_ptr->isotope_unknowns[j].elt_name !=
|
|
string_hsave(kit->second.Get_elt_name().c_str()) ||
|
|
inv_ptr->isotope_unknowns[j].isotope_number !=
|
|
kit->second.Get_isotope_number())
|
|
continue;
|
|
d1 = kit->second.Get_ratio();
|
|
d2 = inv_delta1[col_isotopes +
|
|
i * inv_ptr->isotope_unknowns.size() +
|
|
j] / inv_delta1[i];
|
|
d3 = d1 + d2;
|
|
kit->second.Set_ratio(d3);
|
|
}
|
|
}
|
|
}
|
|
|
|
set_initial_solution(-6, -7);
|
|
temp = pr.all;
|
|
pr.all = FALSE;
|
|
temp_punch = pr.punch;
|
|
pr.punch = FALSE;
|
|
phrq_io->Set_punch_on(false);
|
|
initial_solutions(FALSE);
|
|
pr.all = temp;
|
|
pr.punch = temp_punch;
|
|
phrq_io->Set_punch_on(pr.punch == TRUE);
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, -7);
|
|
|
|
/* Header */
|
|
std::string string;
|
|
if (solution_ptr_orig->Get_description().size() > 0)
|
|
{
|
|
fprintf(netpath_file, "%d. %s\n", count_inverse_models,
|
|
solution_ptr_orig->Get_description().c_str());
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%d. Solution %d\n", count_inverse_models,
|
|
solution_ptr_orig->Get_n_user());
|
|
}
|
|
|
|
/* bookkeeping */
|
|
count_pat_solutions++;
|
|
solnmap[count_current_solutions][0] = solution_ptr_orig->Get_n_user();
|
|
solnmap[count_current_solutions][1] = count_pat_solutions;
|
|
count_current_solutions++;
|
|
|
|
/* Dump info to .pat file */
|
|
print_total_pat(netpath_file, "C", "C");
|
|
print_total_pat(netpath_file, "S", "S");
|
|
print_total_pat(netpath_file, "Ca", "CA");
|
|
print_total_pat(netpath_file, "Al", "AL");
|
|
print_total_pat(netpath_file, "Mg", "MG");
|
|
print_total_pat(netpath_file, "Na", "NA");
|
|
print_total_pat(netpath_file, "K", "K");
|
|
print_total_pat(netpath_file, "Cl", "CL");
|
|
print_total_pat(netpath_file, "F", "F");
|
|
print_total_pat(netpath_file, "Si", "SI");
|
|
print_total_pat(netpath_file, "Br", "BR");
|
|
print_total_pat(netpath_file, "B", "B");
|
|
print_total_pat(netpath_file, "Ba", "BA");
|
|
print_total_pat(netpath_file, "Li", "LI");
|
|
print_total_pat(netpath_file, "Sr", "SR");
|
|
print_total_pat(netpath_file, "Fe", "FE");
|
|
print_total_pat(netpath_file, "Mn", "MN");
|
|
print_total_pat(netpath_file, "N", "N");
|
|
print_total_pat(netpath_file, "P", "P");
|
|
fprintf(netpath_file, "%14g # TEMP\n", (double) solution_ptr->Get_tc());
|
|
print_total_pat(netpath_file, "S(-2)", "H2S");
|
|
print_total_pat(netpath_file, "S(6)", "SO4");
|
|
|
|
/* N15 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
std::map < std::string, cxxSolutionIsotope >::iterator kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "15N") != NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # N15\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # N15\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/* RS of N */
|
|
sum = 0;
|
|
sum =
|
|
total("N(-3)") * -3 + total("N(0)") * 0 + total("N(3)") * 3 +
|
|
total("N(5)") * 5;
|
|
sum1 = total("N(-3)") + total("N(0)") + total("N(3)") + total("N(5)");
|
|
if (sum1 == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # RS of N\n", (double) sum1);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # RS of N\n",
|
|
(double) (sum / sum1));
|
|
}
|
|
|
|
/* DOX */
|
|
print_total_pat(netpath_file, "O(0)", "DOX");
|
|
|
|
/*HCO3 */
|
|
d = 1000 * sum_match_species("*HCO3*", "C");
|
|
if (d == 0.0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # HCO3\n", (double) d);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # HCO3\n", (double) d);
|
|
}
|
|
|
|
/* pH */
|
|
fprintf(netpath_file, "%14g # PH\n", (double) solution_ptr->Get_ph());
|
|
|
|
/*H2CO3* */
|
|
d = 1000 * (molality("H2CO3") + molality("CO2"));
|
|
if (d == 0.0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # H2CO3\n", (double) d);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # H2CO3\n", (double) d);
|
|
}
|
|
|
|
/*CO3 */
|
|
d = sum_match_species("*CO3*", "C");
|
|
d -= sum_match_species("*HCO3*", "C");
|
|
d *= 1000.0;
|
|
if (d == 0.0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # CO3\n", (double) d);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # CO3\n", (double) d);
|
|
}
|
|
|
|
/* CARBONATES */
|
|
print_total_pat(netpath_file, "C(4)", "CARBONATES");
|
|
print_total_pat(netpath_file, "Fe(2)", "FE2+");
|
|
print_total_pat(netpath_file, "Fe(3)", "FE3+");
|
|
print_total_pat(netpath_file, "Mn(2)", "MN2+");
|
|
print_total_pat(netpath_file, "Mn(3)", "MN3+");
|
|
print_total_pat(netpath_file, "Mn(6)", "MN6+");
|
|
print_total_pat(netpath_file, "Mn(7)", "MN7+");
|
|
print_total_pat(netpath_file, "C(-4)", "CH4");
|
|
print_total_pat(netpath_file, "Doc", "DOC");
|
|
|
|
/*RS OF DOC */
|
|
fprintf(netpath_file, "%14g* # RS OF DOC\n", 0.0);
|
|
|
|
/* Blank */
|
|
print_total_pat(netpath_file, "Blank", "BLANK");
|
|
|
|
/*C13 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "13C") != NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # C13\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # C13\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*C14 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "14C") != NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # C14\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # C14\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*SR87 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "87Sr") !=
|
|
NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # SR87\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # SR87\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*D*/ sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "2H") != NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # D\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # D\n", (double) (sum_iso / sum));
|
|
}
|
|
|
|
/*O-18 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "18O") != NULL)
|
|
{
|
|
if (strcmp(kit->second.Get_elt_name().c_str(), "O(-2)") == 0)
|
|
{
|
|
d = solution_ptr->Get_total_o() - total("O(0)");
|
|
}
|
|
else if (strcmp(kit->second.Get_elt_name().c_str(), "H(1)")
|
|
== 0)
|
|
{
|
|
d = solution_ptr->Get_total_h() - total("H(0)");
|
|
}
|
|
else
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
}
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # O-18\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # O-18\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*TRITIUM*/ sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "3H") != NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # TRITIUM\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # TRITIUM\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*34SSO4 */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "34S(6)") !=
|
|
NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # 34SSO4\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # 34SSO4\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/*34SH2S */
|
|
sum_iso = 0;
|
|
sum = 0;
|
|
kit = solution_ptr->Get_isotopes().begin();
|
|
for ( ; kit != solution_ptr->Get_isotopes().end(); kit++)
|
|
{
|
|
if (strstr(kit->second.Get_isotope_name().c_str(), "34S(-2)") !=
|
|
NULL)
|
|
{
|
|
d = total(kit->second.Get_elt_name().c_str());
|
|
sum_iso += kit->second.Get_ratio() * d;
|
|
sum += d;
|
|
}
|
|
}
|
|
if (sum == 0)
|
|
{
|
|
fprintf(netpath_file, "%14g* # 34SH2S\n", (double) sum);
|
|
}
|
|
else
|
|
{
|
|
fprintf(netpath_file, "%14g # 34SH2S\n",
|
|
(double) (sum_iso / sum));
|
|
}
|
|
|
|
/* Well number */
|
|
fprintf(netpath_file, "%14d # Well number\n",
|
|
count_pat_solutions);
|
|
}
|
|
my_array = array_save;
|
|
delta = l_delta_save;
|
|
count_unknowns = count_unknowns_save;
|
|
max_row_count = max_row_count_save;
|
|
max_column_count = max_column_count_save;
|
|
|
|
/*
|
|
* Open model file
|
|
*/
|
|
std::string string;
|
|
string = inv_ptr->pat;
|
|
replace(".pat", "", string);
|
|
trim(string);
|
|
std::string string1 = sformatf("%s-%d.mod", string.c_str(), count_inverse_models);
|
|
model_file = fopen(string1.c_str(), "w");
|
|
if (model_file == NULL)
|
|
{
|
|
error_string = sformatf( "Can`t open file, %s.", string.c_str());
|
|
error_msg(error_string, STOP);
|
|
}
|
|
add_to_file("model.fil", string1.c_str());
|
|
|
|
|
|
/*
|
|
* Write header
|
|
*/
|
|
fprintf(model_file, "%s\n", string.c_str());
|
|
|
|
/*
|
|
* Write well numbers
|
|
*/
|
|
for (i = 0; i < count_current_solutions; i++)
|
|
{
|
|
fprintf(model_file, "%3d", solnmap[i][1]);
|
|
}
|
|
fprintf(model_file, "\n");
|
|
/*
|
|
* Write elements
|
|
*/
|
|
xsolution_zero();
|
|
for (j = 0; j < (int)master.size(); j++)
|
|
{
|
|
master[j]->in = FALSE;
|
|
}
|
|
for (j = 0; j < inv_ptr->elts.size(); j++)
|
|
{
|
|
master_ptr = inv_ptr->elts[j].master;
|
|
master_ptr = master_ptr->elt->primary;
|
|
if (strcmp(master_ptr->elt->name, "Alkalinity") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "H") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "O") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "X") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "E") == 0)
|
|
continue;
|
|
master_ptr->in = TRUE;
|
|
}
|
|
for (j = 0; j < (int)master.size(); j++)
|
|
{
|
|
if (master[j]->in == TRUE)
|
|
{
|
|
string = master[j]->elt->name;
|
|
Utilities::str_toupper(string);
|
|
fprintf(model_file, " %-2s", string.c_str());
|
|
}
|
|
}
|
|
fprintf(model_file, " %-2s", "RS");
|
|
/*
|
|
* Add isotope mole balance
|
|
*/
|
|
for (j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
string = sformatf("%d%s", (int) inv_ptr->isotopes[j].isotope_number,
|
|
inv_ptr->isotopes[j].elt_name);
|
|
if (strcmp(string.c_str(), "13C") == 0)
|
|
fprintf(model_file, " %-2s", "I1");
|
|
if (strcmp(string.c_str(), "14C") == 0)
|
|
fprintf(model_file, " %-2s", "I2");
|
|
if (strcmp(string.c_str(), "34S") == 0)
|
|
fprintf(model_file, " %-2s", "I3");
|
|
if (strcmp(string.c_str(), "87Sr") == 0)
|
|
fprintf(model_file, " %-2s", "I4");
|
|
if (strcmp(string.c_str(), "15N") == 0)
|
|
fprintf(model_file, " %-2s", "I9");
|
|
if (strcmp(string.c_str(), "2H") == 0)
|
|
fprintf(model_file, " %-2s", "D ");
|
|
if (strcmp(string.c_str(), "3H") == 0)
|
|
fprintf(model_file, " %-2s", "TR");
|
|
if (strcmp(string.c_str(), "18O") == 0)
|
|
fprintf(model_file, " %-2s", "18");
|
|
}
|
|
|
|
/* end of element line */
|
|
fprintf(model_file, "\n");
|
|
|
|
/*
|
|
* Write phase information
|
|
*/
|
|
for (size_t i = 0; i < inv_ptr->phases.size(); i++)
|
|
{
|
|
size_t j = col_phases + i;
|
|
/* skip if not in model */
|
|
/* if (equal (inv_delta1[j], 0.0, toler) == TRUE) continue;*/
|
|
|
|
/* Do not include Na exchange phase */
|
|
if (strcmp_nocase(inv_ptr->phases[i].name, "NaX") == 0)
|
|
continue;
|
|
/*
|
|
* Determine if exchange reaction
|
|
*/
|
|
exch = FALSE;
|
|
for (next_elt = &inv_ptr->phases[i].phase->next_elt[0];
|
|
next_elt->elt != NULL; next_elt++)
|
|
{
|
|
if (strcmp(next_elt->elt->name, "X") == 0)
|
|
{
|
|
exch = TRUE;
|
|
break;
|
|
}
|
|
}
|
|
/*
|
|
* Write phase name and constraints
|
|
*/
|
|
string = inv_ptr->phases[i].name;
|
|
string = string.substr(0,8);
|
|
string = Utilities::pad_right(string, 8);
|
|
if (inv_ptr->phases[i].force == TRUE)
|
|
{
|
|
string += 'F';
|
|
}
|
|
else
|
|
{
|
|
string += ' ';
|
|
}
|
|
switch (inv_ptr->phases[i].constraint)
|
|
{
|
|
case EITHER:
|
|
string += ' ';
|
|
break;
|
|
case PRECIPITATE:
|
|
if (exch == TRUE)
|
|
{
|
|
string += '+';
|
|
}
|
|
else
|
|
{
|
|
string += '-';
|
|
}
|
|
break;
|
|
case DISSOLVE:
|
|
if (exch == TRUE)
|
|
{
|
|
string += '-';
|
|
}
|
|
else
|
|
{
|
|
string += '+';
|
|
}
|
|
break;
|
|
}
|
|
fprintf(model_file, "%-10s", string.c_str());
|
|
/*
|
|
* Write stoichiometry
|
|
*/
|
|
for (next_elt = &inv_ptr->phases[i].phase->next_elt[0];
|
|
next_elt->elt != NULL; next_elt++)
|
|
{
|
|
f = 1.0;
|
|
if (exch == TRUE)
|
|
f = -1.0;
|
|
master_ptr = next_elt->elt->primary;
|
|
if (strcmp(master_ptr->elt->name, "Alkalinity") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "H") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "O") == 0)
|
|
continue;
|
|
if (strcmp(master_ptr->elt->name, "E") == 0)
|
|
continue;
|
|
string = master_ptr->elt->name;
|
|
|
|
if (strcmp(master_ptr->elt->name, "X") == 0)
|
|
{
|
|
string = "Na";
|
|
f = 1.0;
|
|
}
|
|
Utilities::str_toupper(string);
|
|
fprintf(model_file, " %-2s%12.7f", string.c_str(),
|
|
(double) (next_elt->coef * f));
|
|
}
|
|
/*
|
|
* Calculate RS
|
|
*/
|
|
std::string token;
|
|
sum = 0;
|
|
for (rxn_ptr = &inv_ptr->phases[i].phase->rxn_s.token[0] + 1;
|
|
rxn_ptr->s != NULL; rxn_ptr++)
|
|
{
|
|
if (rxn_ptr->s == s_hplus)
|
|
continue;
|
|
if (rxn_ptr->s == s_h2o)
|
|
continue;
|
|
if (rxn_ptr->s->secondary == NULL && rxn_ptr->s != s_eminus)
|
|
continue;
|
|
if (rxn_ptr->s == s_o2)
|
|
{
|
|
sum += 4 * rxn_ptr->coef;
|
|
}
|
|
else if (rxn_ptr->s == s_h2)
|
|
{
|
|
sum += -2 * rxn_ptr->coef;
|
|
}
|
|
else if (rxn_ptr->s == s_eminus)
|
|
{
|
|
sum += -1 * rxn_ptr->coef;
|
|
}
|
|
else
|
|
{
|
|
string = rxn_ptr->s->secondary->elt->name;
|
|
replace("(", " ", string);
|
|
replace(")", " ", string);
|
|
std::string::iterator b = string.begin();
|
|
std::string::iterator e = string.end();
|
|
CParser::copy_token(token, b, e);
|
|
CParser::copy_token(string1, b, e);
|
|
(void)sscanf(string1.c_str(), SCANFORMAT, &f);
|
|
sum += f * rxn_ptr->coef;
|
|
}
|
|
}
|
|
if (sum != 0.0)
|
|
fprintf(model_file, " %-2s%12.7f", "RS", (double) sum);
|
|
/*
|
|
* Add isotopes
|
|
*/
|
|
|
|
for (k = 0; k < inv_ptr->phases[i].isotopes.size(); k++)
|
|
{
|
|
std::vector<class isotope>& isotope_ref = inv_ptr->phases[i].isotopes;
|
|
d1 = isotope_ref[k].ratio;
|
|
for (j = 0; j < inv_ptr->isotopes.size(); j++)
|
|
{
|
|
if ((inv_ptr->isotopes[j].elt_name != isotope_ref[k].elt_name)
|
|
|| (inv_ptr->isotopes[j].isotope_number !=
|
|
isotope_ref[k].isotope_number))
|
|
continue;
|
|
break;
|
|
}
|
|
d2 = 0.0;
|
|
if (j < inv_ptr->isotopes.size())
|
|
{
|
|
column = col_phase_isotopes + i * inv_ptr->isotopes.size() + j;
|
|
if (inv_delta1[col_phases + i] != 0.0)
|
|
{
|
|
d2 = inv_delta1[column] / inv_delta1[col_phases + i];
|
|
}
|
|
}
|
|
d3 = d1 + d2;
|
|
string = sformatf("%d%s", (int)isotope_ref[k].isotope_number,
|
|
isotope_ref[k].elt_name);
|
|
if (strcmp(string.c_str(), "13C") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "I1", (double) d3);
|
|
if (strcmp(string.c_str(), "14C") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "I2", (double) d3);
|
|
if (strcmp(string.c_str(), "34S") == 0)
|
|
{
|
|
fprintf(model_file, " %-2s%12.7f", "I3", (double) d3);
|
|
fprintf(model_file, " %-2s%12.7f", "I7", 0.0);
|
|
}
|
|
if (strcmp(string.c_str(), "87Sr") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "I4", (double) d3);
|
|
if (strcmp(string.c_str(), "15N") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "I9", (double) d3);
|
|
if (strcmp(string.c_str(), "2H") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "D ", (double) d3);
|
|
if (strcmp(string.c_str(), "3H") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "TR", (double) d3);
|
|
if (strcmp(string.c_str(), "18O") == 0)
|
|
fprintf(model_file, " %-2s%12.7f", "18", (double) d3);
|
|
}
|
|
/*
|
|
* Done with stoichiometry
|
|
*/
|
|
fprintf(model_file, "\n");
|
|
}
|
|
fprintf(model_file, "\n");
|
|
/*
|
|
* Write extra stuff at bottom
|
|
*/
|
|
/*
|
|
(Iflag(i),i=2,6), (P(i),i=1,2), (Isdocrs(i),i=0,5), Disalong,
|
|
(C14dat(i),i=1,5), (Usera(i),i=1,2),
|
|
(C14dat(i),i=8,9), i10, i11, (C14dat(i),i=12,13),
|
|
(Dbdata(Well(0),i),i=44,47), Dbdata(Well(0),49),
|
|
((Dbdata(Well(iwell),i),i=44,47),Dbdata(Well(iwell),49),Usera(iwell),iwell=1,Iflag(1)+1)
|
|
9030 FORMAT (5(I2),2(F8.4),6(I1),F6.3,/,
|
|
7(F8.3),/,
|
|
2(F8.3),2(F8.0), 2(F8.3),/,
|
|
5(F8.3),/,
|
|
7(6(F8.3),/))
|
|
*/
|
|
/* iflags */
|
|
/*fprintf(model_file,"%2d", i); */ /* not written, 1, mixing, number of mixing wells -1 */
|
|
fprintf(model_file, "%2d", 3); /* 2, exchange */
|
|
i = 0;
|
|
if (inv_ptr->isotopes.size() > 0)
|
|
i = 1;
|
|
fprintf(model_file, "%2d", i); /* 3, Rayleigh */
|
|
fprintf(model_file, "%2d", 1); /* 4, A0 model */
|
|
fprintf(model_file, "%2d", 0); /* 5, Mook/Deines */
|
|
fprintf(model_file, "%2d", 0); /* 6, Evaporation/Dilution */
|
|
/* p */
|
|
fprintf(model_file, "%8.4f%8.4f", 1.0, 0.0); /* P(1),(2) fraction of CO2, fraction of Ca in exch */
|
|
fprintf(model_file, "000000"); /* isdocrs(0,5) doc, rs? */
|
|
fprintf(model_file, "%6.3f\n", 1.0); /* disalong */
|
|
fprintf(model_file,
|
|
" 0.000 100.000 0.000 0.000 -25.000 100.000 100.000\n");
|
|
fprintf(model_file,
|
|
" 0.000 0.000 0. 0. 0.000 0.000 \n");
|
|
/* Final well data */
|
|
fprintf(model_file, " -40.000 -25.000 0.000 0.000 0.000\n");
|
|
/* other wells */
|
|
for (i = 0; i < count_current_solutions - 1; i++)
|
|
{
|
|
fprintf(model_file,
|
|
" -40.000 -25.000 0.000 0.000 0.000 100.000\n");
|
|
}
|
|
fprintf(model_file, "\n");
|
|
fclose(model_file);
|
|
state = INVERSE;
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
void Phreeqc::
|
|
print_total_pat(FILE * l_netpath_file, const char *elt, const char *string)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
LDBLE d;
|
|
d = 1000.0 * total(elt);
|
|
if (strcmp(elt,"O(0)") == 0)
|
|
{
|
|
d = d/2.0;
|
|
}
|
|
if (d == 0)
|
|
{
|
|
fprintf(l_netpath_file, "%14g%1s # %s\n", (double) d, "*", string);
|
|
}
|
|
else
|
|
{
|
|
fprintf(l_netpath_file, "%14g%1s # %s\n", (double) d, " ", string);
|
|
}
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
set_initial_solution(int n_user_old, int n_user_new)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
cxxSolution *solution_ptr;
|
|
Utilities::Rxn_copy(Rxn_solution_map, n_user_old, n_user_new);
|
|
Rxn_new_solution.insert(n_user_new);
|
|
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, n_user_new);
|
|
solution_ptr->Set_new_def(true);
|
|
if (solution_ptr->Get_initial_data() == NULL)
|
|
solution_ptr->Create_initial_data();
|
|
solution_ptr->Set_n_user_end(n_user_new);
|
|
cxxNameDouble::iterator jit = solution_ptr->Get_totals().begin();
|
|
for ( ; jit != solution_ptr->Get_totals().end(); jit++)
|
|
{
|
|
cxxISolutionComp temp_comp;
|
|
temp_comp.Set_description(jit->first.c_str());
|
|
temp_comp.Set_input_conc(jit->second / solution_ptr->Get_mass_water());
|
|
temp_comp.Set_units("Mol/kgw");
|
|
solution_ptr->Get_initial_data()->Get_comps()[jit->first.c_str()] = temp_comp;
|
|
}
|
|
return (OK);
|
|
}
|
|
/* ---------------------------------------------------------------------- */
|
|
int Phreeqc::
|
|
add_to_file(const char *filename, const char *string)
|
|
/* ---------------------------------------------------------------------- */
|
|
{
|
|
FILE *model_file;
|
|
char c;
|
|
int i;
|
|
char string_line[MAX_LINE];
|
|
|
|
model_file = fopen(filename, "r");
|
|
if (model_file == NULL)
|
|
{
|
|
model_file = fopen(filename, "w");
|
|
}
|
|
if (model_file == NULL)
|
|
{
|
|
error_string = sformatf( "Can`t open file, %s.", filename);
|
|
error_msg(error_string, STOP);
|
|
#if !defined(R_SO)
|
|
exit(4);
|
|
#endif
|
|
}
|
|
i = 0;
|
|
/*
|
|
* Read each line of file, check if equal to string; if not, append string to file
|
|
*/
|
|
for (;;)
|
|
{
|
|
c = getc(model_file);
|
|
if (c != EOF && c != '\n' && i != MAX_LINE)
|
|
{
|
|
string_line[i] = c;
|
|
i++;
|
|
continue;
|
|
}
|
|
if (i >= MAX_LINE)
|
|
{
|
|
string_line[MAX_LINE - 1] = '\0';
|
|
error_string = sformatf( "File name in %s is greater than %d characters: %s\n", filename, MAX_LINE, string_line);
|
|
warning_msg(error_string);
|
|
}
|
|
else
|
|
{
|
|
string_line[i] = '\0';
|
|
}
|
|
/* new line or eof */
|
|
string_trim(string_line);
|
|
if (strcmp(string_line, string) == 0)
|
|
{
|
|
fclose(model_file);
|
|
return (OK);
|
|
}
|
|
/* eof, add line */
|
|
if (c == EOF)
|
|
{
|
|
fclose(model_file);
|
|
model_file = fopen(filename, "a");
|
|
if (model_file)
|
|
{
|
|
fprintf(model_file, "%s\n", string);
|
|
fclose(model_file);
|
|
}
|
|
else
|
|
{
|
|
error_string = sformatf("Could not open netpath model file: %s\n", filename);
|
|
error_msg(error_string, STOP);
|
|
}
|
|
return (OK);
|
|
}
|
|
/* new line */
|
|
i = 0;
|
|
}
|
|
}
|
|
|