max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 16:44:49 +01:00
Code Issues Packages Projects Releases Wiki Activity

merged in http://internalbrr.cr.usgs.gov/svn_GW/phreeqcpp/branches/ErrorHandling 6269-6445

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/trunk@6446 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
Scott R Charlton 2012-04-11 02:21:16 +00:00
parent e972c511f0
commit 495a220070
18 changed files with 630 additions and 816 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
ISolution.cxx
View File

@ -18,12 +18,8 @@
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
//static std::map<int, cxxISolution> ss_map;
//std::map<int, cxxISolution>& cxxISolution::s_map = ss_map;
cxxISolution::cxxISolution(PHRQ_io *io)
:
//cxxSolution(io),
units("mMol/kgw")
{
default_pe = "pe";
@ -31,16 +27,8 @@ units("mMol/kgw")
pe_reactions[default_pe] = temp_pe_reactions;
}
//cxxISolution::cxxISolution(const cxxISolution *is_old)
//{
// units = is_old->units;
// comps = is_old->comps;
// pe_reactions = is_old->pe_reactions;
// default_pe = is_old->default_pe;
//}
cxxISolution::~cxxISolution()
{
//// ToDo //pe_data_free(this->pes);
}
#ifdef SKIP_OR_MOVE_TO_STRUCTURES
@ -135,244 +123,6 @@ cxxISolution::ConvertUnits(PHREEQC_PTR_ARG)
}
#endif
#ifdef SKIP
cxxISolution & cxxISolution::read(CParser & parser)
{
static
std::vector <
std::string >
vopts;
if (vopts.empty())
{
vopts.reserve(11);
vopts.push_back("temp"); // 0
vopts.push_back("temperature"); // 1
vopts.push_back("dens"); // 2
vopts.push_back("density"); // 3
vopts.push_back("units"); // 4
vopts.push_back("redox"); // 5
vopts.push_back("ph"); // 6
vopts.push_back("pe"); // 7
vopts.push_back("unit"); // 8
vopts.push_back("isotope"); // 9
vopts.push_back("water"); // 10
}
// const int count_opt_list = vopts.size();
cxxISolution
numkey;
// Read solution number and description
numkey.read_number_description(parser);
// Malloc space for solution data
//// g_solution_map[numkey.n_user()] = numkey;
s_map[numkey.n_user()] = numkey;
std::istream::pos_type ptr;
std::istream::pos_type next_char;
std::string token;
CParser::TOKEN_TYPE j;
cxxISolution & sol = s_map[numkey.n_user()];
int
default_pe = 0;
for (;;)
{
int
opt = parser.get_option(vopts, next_char);
if (opt == CParser::OPTION_DEFAULT)
{
ptr = next_char;
if (parser.copy_token(token, ptr) == CParser::TT_DIGIT)
{
opt = 9;
}
}
switch (opt)
{
case CParser::OPTION_EOF:
break;
case CParser::OPTION_KEYWORD:
break;
case CParser::OPTION_ERROR:
opt = CParser::OPTION_EOF;
CParser::error_msg("Unknown input in SOLUTION keyword.",
PHRQ_io::OT_CONTINUE);
CParser::error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE);
break;
case 0: // temp
case 1: // temperature
if (!(parser.get_iss() >> sol.tc))
{
sol.tc = 25;
}
break;
case 2: // dens
case 3: // density
parser.get_iss() >> sol.density;
break;
case 4: // units
case 8: // unit
if (parser.copy_token(token, next_char) == CParser::TT_EMPTY)
break;
if (parser.check_units(token, false, false, sol.units, true) ==
CParser::OK)
{
sol.units = token;
}
else
{
parser.incr_input_error();
}
break;
case 5: // redox
if (parser.copy_token(token, next_char) == CParser::TT_EMPTY)
break;
if (parser.parse_couple(token) == CParser::OK)
{
default_pe = cxxPe_Data::store(sol.pe, token);
}
else
{
parser.incr_input_error();
}
break;
case 6: // ph
{
cxxISolutionComp
conc;
if (conc.read(parser, sol) == cxxISolutionComp::ERROR)
{
parser.incr_input_error();
break;
}
sol.ph = conc.get_input_conc();
if (conc.get_equation_name().empty())
{
break;
}
conc.Set_description("H(1)");
sol.add(conc);
}
break;
case 7: // pe
{
cxxISolutionComp
conc;
if (conc.read(parser, sol) == cxxISolutionComp::ERROR)
{
parser.incr_input_error();
break;
}
sol.solution_pe = conc.get_input_conc();
if (conc.get_equation_name().empty())
{
break;
}
conc.Set_description("E");
sol.add(conc);
}
break;
case 9: // isotope
{
cxxIsotope
isotope;
if (isotope.read(parser) == cxxIsotope::OK)
{
sol.add(isotope);
}
}
break;
case 10: // water
j = parser.copy_token(token, next_char);
if (j == CParser::TT_EMPTY)
{
sol.mass_water = 1.0;
}
else if (j != CParser::TT_DIGIT)
{
parser.incr_input_error();
parser.
error_msg
("Expected numeric value for mass of water in solution.",
PHRQ_io::OT_CONTINUE);
}
else
{
std::istringstream(token) >> sol.mass_water;
}
break;
case CParser::OPTION_DEFAULT:
{
// Read concentration
cxxISolutionComp
conc;
if (conc.read(parser, sol) == cxxISolutionComp::ERROR)
{
parser.incr_input_error();
}
else
{
sol.add(conc);
}
}
break;
}
if (opt == CParser::OPTION_EOF || opt == CParser::OPTION_KEYWORD)
break;
}
//
// fix up default units and default pe
//
std::string token1;
std::vector < cxxISolutionComp >::iterator iter = sol.totals.begin();
for (; iter != sol.totals.end(); ++iter)
{
token = (*iter).Get_description();
Utilities::str_tolower(token);
if ((*iter).get_units().empty())
{
(*iter).set_units(sol.units);
}
else
{
bool
alk = false;
if (token.find("alk") == 0)
alk = true;
token1 = (*iter).get_units();
if (parser.
check_units(token1, alk, true, sol.get_units(),
true) == CParser::ERROR)
{
parser.incr_input_error();
}
else
{
(*iter).set_units(token1);
}
}
if ((*iter).get_n_pe() < 0)
{
(*iter).set_n_pe(default_pe);
}
}
sol.default_pe = default_pe;
return sol;
}
#endif
#ifdef SKIP
void
cxxISolution::dump_xml(std::ostream & os, unsigned int indent) const const

73
Makefile
View File

@ -32,10 +32,20 @@ Release: class_release
release: class_release
Class_release: class_release
Debug_64: class_debug_64
debug_64: class_debug_64
Class_debug_64: class_debug_64
CLASS_DEBUG_DIR = Class_debug
Release_64: class_release_64
release_64: class_release_64
Class_release_64: class_release_64
CLASS_DEBUG_DIR = Class_debug
CLASS_DIR = Class_release
MAKEFILE = Makefile
CLASS_DEBUG_64_DIR = Class_debug_64
CLASS_64_DIR = Class_release_64
MAKEFILE = Makefile
# -----------------------------------------------------------------------------
# fixes shared object lookup error(SIGFPE floating point exception)
@ -54,6 +64,15 @@ class_debug:
class_release:
mkdir -p $(CLASS_DIR)
cd $(CLASS_DIR); $(MAKE) -r -f ../$(MAKEFILE) CFG=CLASS_RELEASE $(PROGRAM)
.PHONY : Class_debug_64
class_debug_64:
mkdir -p $(CLASS_DEBUG_64_DIR)
cd $(CLASS_DEBUG_64_DIR); $(MAKE) -r -f ../$(MAKEFILE) CFG=CLASS_DEBUG_64 $(PROGRAM)
.PHONY : Class_release_64
class_release_64:
mkdir -p $(CLASS_64_DIR)
cd $(CLASS_64_DIR); $(MAKE) -r -f ../$(MAKEFILE) CFG=CLASS_RELEASE_64 $(PROGRAM)
# Recursive make begins here
#
@ -89,23 +108,23 @@ class_release:
# -----------------------------------------------------------------------------
# #define gmp for inverse modeling
# comment the following lines to remove multiprecision option
INVERSE_CL1MP=TRUE
ifdef INVERSE_CL1MP
DEFINE_INVERSE_CL1MP=-DINVERSE_CL1MP
CL1MP_OBJS=cl1mp.o
# CL1MP_LIB=-lgmp
CL1MP_LIB=/z/parkplace/usr/lib/libgmp.a
endif
# -----------------------------------------------------------------------------
#efence for debugging
EFENCE_LIB=-L$(HOME)/packages/efence
# -----------------------------------------------------------------------------
# 2 Versions
# 4 Versions
# -----------------------------------------------------------------------------
ifeq ($(CFG), CLASS_DEBUG)
DEFINES = -DUSE_PHRQ_ALLOC $(DEFINE_INVERSE_CL1MP) -DPHREEQC2
INVERSE_CL1MP=TRUE
ifdef INVERSE_CL1MP
DEFINE_INVERSE_CL1MP=-DINVERSE_CL1MP
CL1MP_OBJS=cl1mp.o
# CL1MP_LIB=-lgmp
CL1MP_LIB=/z/parkplace/usr/lib/libgmp.a
endif
DEFINES = -DUSE_PHRQ_ALLOC $(DEFINE_INVERSE_CL1MP) # -DPHREEQC2
VPATH = ..:../phreeqc
INCLUDES = -I../phreeqc -I..
CXX = g++
@ -115,12 +134,38 @@ ifeq ($(CFG), CLASS_DEBUG)
endif
ifeq ($(CFG), CLASS_RELEASE)
DEFINES = -DNDEBUG $(DEFINE_INVERSE_CL1MP) -DPHREEQC2
INVERSE_CL1MP=TRUE
ifdef INVERSE_CL1MP
DEFINE_INVERSE_CL1MP=-DINVERSE_CL1MP
CL1MP_OBJS=cl1mp.o
# CL1MP_LIB=-lgmp
CL1MP_LIB=/z/parkplace/usr/lib/libgmp.a
endif
DEFINES = -DNDEBUG $(DEFINE_INVERSE_CL1MP) # -DPHREEQC2
VPATH = ..:../phreeqc
INCLUDES = -I../phreeqc -I..
CXX = g++
CXXFLAGS = -Wall -pedantic -O3 $(DEFINES) $(INCLUDES)
OBJECT_FILES = $(CLASS_FILES) $(COMMON_COBJS) $(COMMON_CXXOBJS) $(CL1MP_OBJS)
LD_FLAGS = -lm ${CL1MP_LIB} ${HASH_STYLE}
endif
ifeq ($(CFG), CLASS_DEBUG_64)
DEFINES = -DUSE_PHRQ_ALLOC # $(DEFINE_INVERSE_CL1MP) # -DPHREEQC2
VPATH = ..:../phreeqc
INCLUDES = -I../phreeqc -I..
CXX = g++
CXXFLAGS = -Wall -g $(DEFINES) $(INCLUDES)
OBJECT_FILES = $(CLASS_FILES) $(COMMON_COBJS) $(COMMON_CXXOBJS) $(CL1MP_OBJS)
LD_FLAGS = -lm ${CL1MP_LIB} ${HASH_STYLE}
endif
ifeq ($(CFG), CLASS_RELEASE_64)
DEFINES = -DNDEBUG # $(DEFINE_INVERSE_CL1MP) # -DPHREEQC2
VPATH = ..:../phreeqc
INCLUDES = -I../phreeqc -I..
CXX = g++
CXXFLAGS = -Wall -pedantic -O3 $(DEFINES) $(INCLUDES)
# CXXFLAGS = -Wall -pedantic -p $(DEFINES) $(INCLUDES)
OBJECT_FILES = $(CLASS_FILES) $(COMMON_COBJS) $(COMMON_CXXOBJS) $(CL1MP_OBJS)
LD_FLAGS = -lm ${CL1MP_LIB} ${HASH_STYLE}
endif
@ -927,7 +972,7 @@ utilities.o: ../phreeqc/utilities.cpp ../Utils.h ../Phreeqc.h \
../ISolutionComp.h
# -----------------------------------------------------------------------------
clean:
rm -rf Class_release Class_debug
rm -rf Class_release Class_debug Class_release_64 Class_debug_64
dependencies:
mkdir -p $(CLASS_DEBUG_DIR)

21
PHRQ_io.cpp
View File

@ -13,11 +13,11 @@
PHRQ_io::
PHRQ_io(void)
{
output_ostream = NULL;
log_ostream = NULL;
output_ostream = NULL;
log_ostream = NULL;
punch_ostream = NULL;
#ifdef ERROR_OSTREAM
error_ostream = NULL;
error_ostream = NULL;
#else
error_file = NULL;
#endif
@ -237,7 +237,7 @@ warning_msg(const char *err_str)
err_stdstr.append("\n");
screen_msg(err_stdstr.c_str());
error_ostream->flush();
}
}
std::ostringstream warn_str;
warn_str << err_str << "\n";
log_msg(warn_str.str().c_str());
@ -365,7 +365,7 @@ warning_msg(const char *err_str)
err_stdstr.append("\n");
screen_msg(err_stdstr.c_str());
error_flush();
}
}
std::ostringstream warn_str;
warn_str << err_str << "\n";
log_msg(warn_str.str().c_str());
@ -709,7 +709,7 @@ get_line(void)
this->accumulated.append("\n");
}
//
// New line character encountered
// New line character encountered
//
return_value = (empty ? LT_EMPTY : LT_OK);
@ -717,7 +717,7 @@ get_line(void)
if (continue_loop) continue;
//
// Determine return_value
//
//
if (return_value == LT_OK)
{
if (check_key(m_line.begin(), m_line.end()))
@ -743,25 +743,26 @@ get_line(void)
std::string::iterator end = m_line.end();
CParser::copy_token(stdtoken, beg, end);
Utilities::str_tolower(stdtoken);
if ((strstr(stdtoken.c_str(),"include$") == stdtoken.c_str()) ||
if ((strstr(stdtoken.c_str(),"include$") == stdtoken.c_str()) ||
(strstr(stdtoken.c_str(),"include_file") == stdtoken.c_str()))
{
std::string file_name;
file_name.assign(beg, end);
file_name = trim(file_name);
if (file_name.size() > 0)
{
std::ifstream *next_stream = new std::ifstream(file_name.c_str(), std::ios_base::in);
if (!next_stream->is_open())
{
std::ostringstream errstr;
errstr << "Could not open include file " << file_name;
errstr << "\n*********** Could not open include file " << file_name << ". ***********\n\n";
delete next_stream;
#if defined(PHREEQCI_GUI)
warning_msg(errstr.str().c_str());
continue;
#else
output_msg(errstr.str().c_str());
error_msg(errstr.str().c_str(), OT_STOP);
#endif
}

21
Phreeqc.h
View File

@ -1099,17 +1099,23 @@ protected:
protected:
PHRQ_io *phrq_io;
PHRQ_io ioInstance;
int same_model;
//int same_temperature;
//int same_pressure;
//bool same_mu;
LDBLE current_tc;
LDBLE current_pa;
LDBLE current_mu;
bool mu_terms_in_logk;
/* ----------------------------------------------------------------------
* STRUCTURES
* ---------------------------------------------------------------------- */
struct model last_model;
int same_model;
int same_temperature;
int same_pressure;
bool same_mu;
struct punch punch;
/* ----------------------------------------------------------------------
* Temperatures
* ---------------------------------------------------------------------- */
@ -1252,6 +1258,7 @@ protected:
LDBLE tc_x;
LDBLE tk_x;
LDBLE patm_x;
LDBLE last_patm_x;
bool numerical_fixed_volume;
bool force_numerical_fixed_volume;
bool switch_numerical;
@ -1573,8 +1580,8 @@ protected:
/* ----------------------------------------------------------------------
* ISOTOPES
* ---------------------------------------------------------------------- */
struct name_coef match_tokens[50];
int count_match_tokens;
//struct name_coef match_tokens[50];
//int count_match_tokens;
int count_master_isotope;
struct master_isotope **master_isotope;
int max_master_isotope;
@ -1614,7 +1621,7 @@ protected:
LDBLE AA_basic, BB_basic, CC, I_m, rho_0;
LDBLE eps_r; // relative dielectric permittivity
#else
LDBLE V_solutes, rho_0, kappa_0, p_sat, ah2o_x0;
LDBLE V_solutes, rho_0, kappa_0, p_sat/*, ah2o_x0*/;
LDBLE eps_r; // relative dielectric permittivity
LDBLE DH_A, DH_B, DH_Av; // Debye-Hueckel A, B and Av
#endif

339
phreeqc/basicsubs.cpp
View File

@ -1423,6 +1423,242 @@ int Phreeqc::
match_elts_in_species(const char *name, const char *mytemplate)
/* ---------------------------------------------------------------------- */
{
/*
* Makes a list of elements with their coefficients, stores elements
* in elt_list at position count_elts. Global variable count_elts is
* updated with each stored element. Also uses static global variable
* paren_count.
*
* Arguments:
* **t_ptr input, point in token string to start looking
* output, is next position to start looking
* coef input, coefficient to multiply subscripts by
*/
int i, i1, l, case_no, match;
char c, c1;
char *ptr, *ptr1;
LDBLE d;
char element[MAX_LENGTH];
char token[MAX_LENGTH], equal_list[MAX_LENGTH];
char token1[MAX_LENGTH], template1[MAX_LENGTH], equal_list1[MAX_LENGTH];
char str[2];
strcpy(token, name);
squeeze_white(token);
replace("(+", "(", token);
if (strstr("token", "++") != NULL)
{
replace("++++++", "+6", token);
replace("+++++", "+5", token);
replace("++++", "+4", token);
replace("+++", "+3", token);
replace("++", "+2", token);
}
if (strstr("token", "--") != NULL)
{
replace("------", "-6", token);
replace("-----", "-5", token);
replace("----", "-4", token);
replace("---", "-3", token);
replace("--", "-2", token);
}
ptr = token;
std::vector< std::pair<std::string, LDBLE> > match_vector;
while ((c = *ptr) != '\0')
{
c1 = *(ptr + 1);
str[0] = c;
str[1] = '\0';
/*
* New element
*/
if (isupper((int) c) || (c == 'e' && c1 == '-') || (c == '['))
{
/*
* Get new element and subscript
*/
if (get_elt(&ptr, element, &l) == ERROR)
{
return (ERROR);
}
if (get_num(&ptr, &d) == ERROR)
{
return (ERROR);
}
std::pair<std::string, LDBLE> pr(element, d);
match_vector.push_back(pr);
}
else
{
std::pair<std::string, LDBLE> pr(str, 1.0);
match_vector.push_back(pr);
ptr += 1;
}
}
/*
* Replace elements with first of equivalent elements
*/
strcpy(template1, mytemplate);
squeeze_white(template1);
ptr = template1;
while (extract_bracket(&ptr, equal_list) == TRUE)
{
replace("{", "", equal_list);
replace("}", "", equal_list);
while (replace(",", " ", equal_list) == TRUE);
ptr1 = equal_list;
/*
* Get first name in a list from template
*/
std::string elt_name;
if (copy_token(elt_name, &ptr1) == EMPTY)
{
error_string = sformatf(
"Expecting a nonempty list of element names in isotope sum. %s",
mytemplate);
error_msg(error_string, CONTINUE);
return (ERROR);
}
std::string replace_name = elt_name;
/*
* Replace in species all equivalent names from template
*/
while (copy_token(elt_name, &ptr1) != EMPTY)
{
for (i = 0; i < (int) match_vector.size(); i++)
{
if (elt_name == match_vector[i].first)
{
match_vector[i].first = replace_name;
}
}
}
}
/*
* Combine contiguous elements
*/
i1 = 0;
for (i = 1; i < (int) match_vector.size(); i++)
{
if ((isupper((int) (match_vector[i].first[0])) != FALSE)
&& (match_vector[i].first == match_vector[i1].first))
{
match_vector[i1].second += match_vector[i].second;
}
else
{
i1++;
match_vector[i1].first = match_vector[i].first;
match_vector[i1].second = match_vector[i].second;
}
}
int count_match_tokens = i1 + 1;
/*
* write out string
*/
token[0] = '\0';
for (i = 0; i < count_match_tokens; i++)
{
strcat(token, match_vector[i].first.c_str());
if (match_vector[i].second != 1.0)
{
sprintf(token1, "%g", (double) match_vector[i].second);
strcat(token, token1);
}
}
/*
* Write a template name using first of equivalent elements
*/
strcpy(template1, mytemplate);
squeeze_white(template1);
ptr = template1;
while (extract_bracket(&ptr, equal_list) == TRUE)
{
strcpy(equal_list1, equal_list);
replace("{", "", equal_list);
replace("}", "", equal_list);
while (replace(",", " ", equal_list) == TRUE);
ptr1 = equal_list;
/*
* Get first name in a list
*/
std::string elt_name;
if (copy_token(elt_name, &ptr1) == EMPTY)
{
error_string = sformatf(
"Expecting a nonempty list of element names in isotope sum. %s",
mytemplate);
error_msg(error_string, CONTINUE);
return (ERROR);
}
replace(equal_list1, elt_name.c_str(), template1);
squeeze_white(template1);
ptr = template1;
}
/*
* Compare string
*/
/* Cases: 0 exact match
* 1 leading wild card
* 2 trailing wild card
* 3 leading and trailing wild card
*/
case_no = 0;
if (template1[0] == '*')
case_no = 1;
l = (int) strlen(template1);
if (template1[l - 1] == '*')
{
if (case_no != 1)
{
case_no = 2;
}
else
{
case_no = 3;
}
}
while (replace("*", "", template1));
match = FALSE;
switch (case_no)
{
case 0:
/* exact match */
if (strcmp(token, template1) == 0)
match = TRUE;
break;
case 1:
/* leading wild card */
if ((ptr = strstr(token, template1)) == NULL)
{
match = FALSE;
}
else
{
if (strcmp(ptr, template1) == 0)
match = TRUE;
}
break;
case 2:
/* trailing wild card */
if (strstr(token, template1) == token)
match = TRUE;
break;
case 3:
/* trailing wild card */
if (strstr(token, template1) != NULL)
match = TRUE;
break;
}
return (match);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
match_elts_in_species(const char *name, const char *mytemplate)
/* ---------------------------------------------------------------------- */
{
/*
* Makes a list of elements with their coefficients, stores elements
* in elt_list at position count_elts. Global variable count_elts is
@ -1648,7 +1884,7 @@ match_elts_in_species(const char *name, const char *mytemplate)
}
return (match);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
extract_bracket(char **string, char *bracket_string)
@ -1673,7 +1909,106 @@ extract_bracket(char **string, char *bracket_string)
*string += 1;
return (TRUE);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
surf_total(const char *total_name, const char *surface_name)
/* ---------------------------------------------------------------------- */
{
/*
* Provides total moles in LDBLE layer
*/
int j;
if (use.Get_surface_ptr() == NULL)
return (0);
/*
* Find surface...
*/
for (j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE)
continue;
std::string token;
token = x[j]->master[0]->elt->name;
replace("_", " ", token);
std::string::iterator b = token.begin();
std::string::iterator e = token.end();
std::string name;
CParser::copy_token(name, b, e);
if (surface_name != NULL)
{
if (strcmp(name.c_str(), surface_name) == 0)
break;
}
else
{
break;
}
}
if (j >= count_unknowns)
return (0);
/*
* find total moles for redox state
*/
LDBLE t = 0;
bool redox = false;
if (strstr(total_name, "(") != NULL)
{
redox = true;
}
for (j = 0; j < count_s_x; j++)
{
if (s_x[j]->type != SURF)
continue;
std::string token;
struct rxn_token *rxn_ptr;
for (rxn_ptr = s_x[j]->rxn_s->token + 1; rxn_ptr->s != NULL; rxn_ptr++)
{
if (redox && rxn_ptr->s->secondary)
{
token = rxn_ptr->s->secondary->elt->name;
}
else if (!redox && rxn_ptr->s->secondary)
{
token = rxn_ptr->s->secondary->elt->primary->elt->name;
}
else if (!redox && rxn_ptr->s->primary)
{
token = rxn_ptr->s->primary->elt->name;
}
else
{
continue;
}
if (strcmp(token.c_str(), total_name) == 0)
{
t += rxn_ptr->coef * s_x[j]->moles;
break;
}
else
// sum all sites in case total_name is a surface name without underscore
{
if (rxn_ptr->s->type == SURF)
{
if (token.find("_") != std::string::npos)
{
token = token.substr(0, token.find("_"));
}
if (strcmp(token.c_str(), total_name) == 0)
{
t += rxn_ptr->coef * s_x[j]->moles;
break;
}
}
}
}
}
return t;
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
surf_total(const char *total_name, const char *surface_name)
@ -1760,7 +2095,7 @@ surf_total(const char *total_name, const char *surface_name)
}
return (0);
}
#endif
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
total(const char *total_name)

66
phreeqc/cl1.cpp
View File

@ -2,7 +2,10 @@
#include <string.h>
#include <math.h>
#include <stdlib.h>
//#define USE_GMP
#if defined(USE_GMP)
#include <gmp.h>
#endif
#include "Phreeqc.h"
#include "phqalloc.h"
@ -31,8 +34,8 @@ cl1(int k, int l, int m, int n,
LDBLE xmin, xmax;
int iout = 0;
// static i runs faster on windows
int i, j;
LDBLE l_z;
register int i, j;
register LDBLE l_z;
int maxit, n1, n2;
LDBLE pivot;
int ia, ii, kk, nk, js;
@ -42,7 +45,20 @@ cl1(int k, int l, int m, int n,
LDBLE zu, zv;
LDBLE tpivot;
int klm, jmn, nkl, jpn;
LDBLE cuv, sum;
LDBLE cuv;
#if defined(USE_GMP)
mpf_set_default_prec(256);
mpf_t sum;
mpf_t z_mpf;
mpf_t q_mpf;
mpf_t d_mpf;
mpf_init(sum);
mpf_init(z_mpf);
mpf_init(q_mpf);
mpf_init(d_mpf);
#else
long double sum;
#endif
int klm1;
int q_dim, cu_dim;
int kode_arg;
@ -298,7 +314,11 @@ cl1(int k, int l, int m, int n,
#endif
for (j = js; j < n1; ++j)
{
sum = 0.;
#if defined(USE_GMP)
mpf_set_d(sum, 0.0);
#else
sum = 0.;
#endif
for (i = 0; i < klm; ++i)
{
ii = q2[i * q_dim + n1].ival;
@ -310,9 +330,21 @@ cl1(int k, int l, int m, int n,
{
l_z = l_cu[ii - 1];
}
sum += q2[i * q_dim + j].dval * l_z;
#if defined(USE_GMP)
mpf_set_d(z_mpf, l_z);
mpf_set_d(q_mpf, q2[i * q_dim + j].dval);
mpf_mul(d_mpf, z_mpf, q_mpf);
mpf_add(sum, sum, d_mpf);
#else
sum += (long double) q2[i * q_dim + j].dval * (long double) l_z;
#endif
}
#if defined(USE_GMP)
q2[klm * q_dim + j].dval = mpf_get_d(sum);
#else
q2[klm * q_dim + j].dval = sum;
#endif
}
for (j = js; j < n; ++j)
{
@ -686,7 +718,11 @@ cl1(int k, int l, int m, int n,
#ifdef DEBUG_CL1
output_msg(sformatf( "L590\n"));
#endif
#if defined(USE_GMP)
mpf_set_d(sum, 0.0);
#else
sum = 0.;
#endif
for (j = 0; j < n; ++j)
{
l_x[j] = 0.;
@ -718,7 +754,13 @@ cl1(int k, int l, int m, int n,
if (ii >= n1 && ii <= nk)
{
/* * DBLE(Q(I,N1)) */
sum += q2[i * q_dim + n].dval;
#if defined(USE_GMP)
mpf_set_d(q_mpf, q2[i * q_dim + n].dval);
mpf_add(sum, sum, d_mpf);
#else
sum += (long double) q2[i * q_dim + n].dval;
#endif
}
}
}
@ -726,7 +768,17 @@ cl1(int k, int l, int m, int n,
#ifdef DEBUG_CL1
output_msg(sformatf( "L640\n"));
#endif
#if defined(USE_GMP)
*l_error = mpf_get_d(sum);
#else
*l_error = sum;
#endif
#if defined(USE_GMP)
mpf_clear(sum);
mpf_clear(z_mpf);
mpf_clear(q_mpf);
mpf_clear(d_mpf);
#endif
/*
* Check calculation
*/

1
phreeqc/cl1mp.cpp
View File

@ -22,6 +22,7 @@ cl1mp(int k, int l, int m, int n,
int *iter, LDBLE * x_arg, LDBLE * res_arg, LDBLE * error_arg,
LDBLE * cu_arg, int *iu, int *s, int check, LDBLE censor_arg)
{
mpf_set_default_prec(256);
/* System generated locals */
union double_or_int
{

25
phreeqc/gases.cpp
View File

@ -41,7 +41,7 @@ setup_fixed_volume_gas(void)
{
gas_unknown = gas_unknowns[0];
}
same_pressure = FALSE;
//same_pressure = FALSE;
return (OK);
}
@ -367,7 +367,8 @@ calc_PR(void)
m_sum += gas_unknowns[i]->moles;
phase_ptr = gas_unknowns[i]->phase;
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
continue;
error_msg("Cannot calculate a mixture of ideal and Peng_Robinson gases,\n please define Tc and Pc for the active gases in PHASES.", STOP);
//continue;
if (!phase_ptr->pr_a)
{
T_c = phase_ptr->t_c;
@ -404,16 +405,16 @@ calc_PR(void)
{
a_aa_sum2 = 0.0;
phase_ptr = gas_unknowns[i]->phase;
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
continue;
//if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
// continue;
b_sum += phase_ptr->fraction_x * phase_ptr->pr_b;
size_t i1;
struct phase *phase_ptr1;
for (i1 = 0; i1 < gas_unknowns.size(); i1++)
{
phase_ptr1 = gas_unknowns[i1]->phase;
if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0)
continue;
//if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0)
// continue;
if (phase_ptr1->fraction_x == 0)
continue;
a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha *
@ -503,7 +504,7 @@ calc_PR(void)
// accept a (possible) whobble in the curve...
// error_msg("No convergence when calculating P in Peng-Robinson.", STOP);
}
if (V_m < v1)
if (V_m < v1 && it < 40)
P = R_TK / (v1 - b_sum) - a_aa_sum / (v1 * (v1 + 2 * b_sum) - b2);
}
}
@ -558,11 +559,11 @@ calc_PR(void)
}
phase_ptr->pr_p = phase_ptr->fraction_x * P;
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
{
phase_ptr->pr_phi = 1;
continue;
}
//if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
//{
// phase_ptr->pr_phi = 1;
// continue;
//}
rz = P * V_m / R_TK;
A = a_aa_sum * P / (R_TK * R_TK);
B = b_sum * P / R_TK;

2
phreeqc/kinetics.cpp
View File

@ -1497,7 +1497,7 @@ set_and_run(int i, int use_mix, int use_kinetics, int nsaver,
else
{
prep();
same_pressure = false;
//same_pressure = false;
k_temp(use.Get_solution_ptr()->Get_tc(), use.Get_solution_ptr()->Get_patm());
set(FALSE);
converge = model();

21
phreeqc/mainsubs.cpp
View File

@ -551,8 +551,13 @@ initialize(void)
print_density = 0;
same_model = FALSE;
same_temperature = FALSE;
same_pressure = FALSE;
//same_temperature = FALSE;
//same_pressure = FALSE;
current_tc = NAN;
current_pa = NAN;
current_mu = NAN;
mu_terms_in_logk = true;
g_iterations = -1;
G_TOL = 1e-8;
save_init(-1);
@ -674,12 +679,12 @@ initialize(void)
pe_step_size_now = pe_step_size;
count_total_steps = 0;
remove_unstable_phases = FALSE;
for (i = 0; i < 50; i++)
{
match_tokens[i].name = NULL;
match_tokens[i].coef = 0;
}
count_match_tokens = 0;
//for (i = 0; i < 50; i++)
//{
// match_tokens[i].name = NULL;
// match_tokens[i].coef = 0;
//}
//count_match_tokens = 0;
initial_solution_isotopes = FALSE;
full_pitzer = FALSE;
always_full_pitzer = FALSE;

111
phreeqc/model.cpp
View File

@ -44,7 +44,7 @@ model(void)
int mass_water_switch_save;
LDBLE old_mu;
same_mu = true;
//same_mu = true;
set_inert_moles();
/* debug_model = TRUE; */
/* debug_prep = TRUE; */
@ -159,13 +159,6 @@ model(void)
}
reset();
}
#ifndef PHREEQC2
if (fabs(old_mu - mu_x) > 1e-3 * mu_x)
{
same_mu = false;
old_mu = mu_x;
}
#endif
gammas(mu_x);
if (molalities(FALSE) == ERROR)
{
@ -429,7 +422,7 @@ check_residuals(void)
}
else
{
if ((residual[i] >= epsilon
if ((residual[i] >= epsilon * 100
&& x[i]->moles > 0.0) /* || stop_program == TRUE */ )
{
remove_unstable_phases = TRUE;
@ -610,9 +603,7 @@ gammas(LDBLE mu)
b = 50.2905 * s3 / c1;
#else
// a and b are calc'd in calc_dielectrics(tc_x, patm_x);
if (!same_mu)
k_temp(tc_x, patm_x);
same_mu = true;
k_temp(tc_x, patm_x);
a = DH_A;
b = DH_B;
#endif
@ -2140,7 +2131,7 @@ mb_gases(void)
gas_unknown->moles > MIN_TOTAL)
{
gas_in = TRUE;
patm_x = gas_phase_ptr->Get_total_p();
//patm_x = gas_phase_ptr->Get_total_p();
}
}
@ -2359,7 +2350,6 @@ molalities(int allow_overflow)
else if (s_x[i]->type == SURF)
{
s_x[i]->moles = exp(s_x[i]->lm * LOG_10); /* formerly * mass water */
}
else
{
@ -2602,16 +2592,6 @@ calc_gas_pressures(void)
V_m = 1.0;
calc_PR(phase_ptrs, 0, tk_x, V_m);
pr_done = true;
if (fabs(gas_phase_ptr->Get_total_p() - patm_x) > 0.01)
{
same_pressure = FALSE;
if (V_m < 0.07)
patm_x = (1. * patm_x + gas_phase_ptr->Get_total_p()) / 2;
else
patm_x = gas_phase_ptr->Get_total_p();
k_temp(tc_x, patm_x);
}
} else
{
gas_phase_ptr->Set_total_p(0);
@ -2696,15 +2676,6 @@ calc_gas_pressures(void)
}
gas_phase_ptr->Set_total_p(1500.0);
}
if (iterations > 1 && fabs(gas_phase_ptr->Get_total_p() - patm_x) > 0.01)
{
same_pressure = FALSE;
if (gas_phase_ptr->Get_total_p() > 1e3)
patm_x = (3 * patm_x + gas_phase_ptr->Get_total_p()) / 4;
else
patm_x = (1 * patm_x + gas_phase_ptr->Get_total_p()) / 2;
k_temp(tc_x, patm_x);
}
}
return (OK);
@ -3495,11 +3466,6 @@ reset(void)
else
{
mu_x += delta[i];
if (patm_x > 1 && fabs(delta[i]) > 1e-3 * mu_x)
{
same_model = false;
k_temp(tc_x, patm_x);
}
}
if (mu_x <= 1e-8)
{
@ -3522,6 +3488,7 @@ reset(void)
(double) delta[i], "delta/c", (double) d));
}
s_h2o->la += d;
ah2o_x = exp(s_h2o->la * LOG_10);
/* pe */
}
else if (x[i]->type == MH)
@ -3633,15 +3600,27 @@ reset(void)
x[i]->moles += delta[i];
if (x[i]->moles < MIN_TOTAL)
x[i]->moles = MIN_TOTAL;
if (x[i] == gas_unknown && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME &&
(gas_phase_ptr->Get_pr_in() || force_numerical_fixed_volume)
&& numerical_fixed_volume && !calculating_deriv)
{
patm_x = gas_phase_ptr->Get_total_p();
k_temp(tc_x, patm_x);
!calculating_deriv)
{
if (debug_model == TRUE)
{
output_msg(sformatf(
"%-10.10s %-9s%10.2e %-9s%10.2e %-6s%10.2e\n",
"Pressure", "old P",
(double) last_patm_x, "new P",
(double) gas_phase_ptr->Get_total_p(), "iter P",
(double) patm_x));
}
patm_x = gas_phase_ptr->Get_total_p();
if (patm_x < 1e-10 && patm_x < p_sat)
{
patm_x = ( 1 * patm_x + p_sat) / 2.0;
}
if (patm_x > 1500)
patm_x = 1500;
}
last_patm_x = patm_x;
}
else if (x[i]->type == SS_MOLES)
{
@ -3844,15 +3823,6 @@ residuals(void)
"Failed Residual %d: %s %d %e\n", iterations,
x[i]->description, i, residual[i]));
converge = FALSE;
#ifndef PHREEQC2
ah2o_x = exp(s_h2o->la * LOG_10);
/* recalculate k's if pressure changes > 0.005 atm by ah2o change... */
if ((use.Get_solution_ptr()->Get_patm() < p_sat && fabs(patm_x - p_sat) > 5e-3) ||
fabs(ah2o_x - ah2o_x0) > (5e-3 / p_sat))
{
k_temp(tc_x, patm_x);
}
#endif
}
}
else if (x[i]->type == MH
@ -3991,6 +3961,15 @@ residuals(void)
x[i]->description, i, residual[i]));
converge = FALSE;
}
if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME &&
(fabs(last_patm_x - patm_x) > 0.001 || fabs(last_patm_x - gas_phase_ptr->Get_total_p()) > 0.001)
&& !calculating_deriv)
{
if (print_fail)
output_msg(sformatf("Failed pressure test %d: %e %e %e\n", iterations, last_patm_x, patm_x,
gas_phase_ptr->Get_total_p()));
converge = FALSE;
}
}
else if (x[i]->type == SS_MOLES)
{
@ -4447,7 +4426,7 @@ set(int initial)
#ifdef PHREEQC2
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
#else
//patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
#endif
/*
@ -4552,9 +4531,6 @@ revise_guesses(void)
LDBLE d;
max_iter = 10;
#ifndef PHREEQC2
same_mu = false;
#endif
gammas(mu_x);
l_iter = 0;
repeat = TRUE;
@ -4711,9 +4687,6 @@ revise_guesses(void)
{
mu_x = 1e-8;
}
#ifndef PHREEQC2
same_mu = false;
#endif
gammas(mu_x);
return (OK);
}
@ -4739,6 +4712,7 @@ sum_species(void)
*/
ph_x = -s_hplus->la;
solution_pe_x = -s_eminus->la;
ah2o_x = exp(s_h2o->la * LOG_10);
#ifdef PHREEQC2
ah2o_x = exp(s_h2o->la * LOG_10);
#else
@ -4932,9 +4906,6 @@ surface_model(void)
{
g_iterations++;
prev_aq_x = mass_water_aq_x;
#ifndef PHREEQC2
same_mu = false;
#endif
k_temp(tc_x, patm_x);
gammas(mu_x);
molalities(TRUE);
@ -4966,7 +4937,7 @@ surface_model(void)
debug_diffuse_layer = TRUE;
}
#ifndef PHREEQC2
same_mu = false;
//same_mu = false;
#else
k_temp(tc_x, patm_x);
#endif
@ -5183,7 +5154,7 @@ numerical_jacobian(void)
calculating_deriv = TRUE;
#ifndef PHREEQC2
same_mu = false;
//same_mu = false;
#endif
gammas(mu_x);
molalities(TRUE);
@ -5249,10 +5220,6 @@ numerical_jacobian(void)
case MU:
d2 = d * mu_x;
mu_x += d2;
#ifndef PHREEQC2
if (fabs(d2 / mu_x) > 1e-3)
same_mu = false;
#endif
gammas(mu_x);
break;
case PP:
@ -5340,10 +5307,6 @@ numerical_jacobian(void)
break;
case MU:
mu_x -= d2;
#ifndef PHREEQC2
if (fabs(d2 / mu_x) > 1e-3)
same_mu = false;
#endif
gammas(mu_x);
break;
case PP:

15
phreeqc/pitzer.cpp
View File

@ -1383,7 +1383,7 @@ set_pz(int initial)
#ifdef PHREEQC2
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
#else
//patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
#endif
/*
* H+, e-, H2O
@ -1645,6 +1645,7 @@ jacobian_pz(void)
Restart:
int pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x);
if (full_pitzer == TRUE)
{
molalities(TRUE);
@ -1682,6 +1683,8 @@ Restart:
d2 = d1;
break;
case PITZER_GAMMA:
if (!full_pitzer)
continue;
x[i]->s->lg += d;
d2 = d;
break;
@ -1710,8 +1713,10 @@ Restart:
x[i]->moles += d2;
break;
case MU:
//continue;
d2 = d * mu_x;
mu_x += d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
break;
case PP:
@ -1765,6 +1770,7 @@ Restart:
break;
case MU:
mu_x -= d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
break;
case GAS_MOLES:
@ -2031,13 +2037,9 @@ check_gammas_pz(void)
converge = FALSE;
}
}
if (fabs(old_mu - mu_x) > 1e-3 * mu_x)
{
same_model = false;
k_temp(tc_x, patm_x);
}
if (fabs(old_mu - mu_x) > tol)
converge = FALSE;
if ((pow((LDBLE) 10.0, s_h2o->la) - AW) > tol)
converge = FALSE;
return converge;
@ -2054,6 +2056,7 @@ gammas_pz()
int i, j;
LDBLE coef;
/* Initialize */
k_temp(tc_x, patm_x);
/*
* Calculate activity coefficients
*/

84
phreeqc/prep.cpp
View File

@ -232,7 +232,6 @@ quick_setup(void)
gas_unknowns[i]->phase->pr_phi = 1.0;
gas_unknowns[i]->phase->pr_p = 0;
}
same_pressure = FALSE;
}
else
{
@ -2772,7 +2771,6 @@ reprep(void)
* Build model again
*/
build_model();
same_model = FALSE;
k_temp(tc_x, patm_x);
return (OK);
@ -3973,7 +3971,8 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
m_sum += phase_ptr->moles_x;
}
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
continue;
error_msg("Cannot calculate a mixture of ideal and Peng_Robinson gases,\n please define Tc and Pc for the active gases in PHASES.", STOP);
//continue;
if (!phase_ptr->pr_a)
{
T_c = phase_ptr->t_c;
@ -4014,14 +4013,14 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
{
a_aa_sum2 = 0.0;
phase_ptr = phase_ptrs[i];
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
continue;
//if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
// continue;
b_sum += phase_ptr->fraction_x * phase_ptr->pr_b;
for (i1 = 0; i1 < n_g; i1++)
{
phase_ptr1 = phase_ptrs[i1];
if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0)
continue;
//if (phase_ptr1->t_c == 0.0 || phase_ptr1->p_c == 0.0)
// continue;
if (phase_ptr1->fraction_x == 0)
continue;
a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha *
@ -4072,7 +4071,7 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
27. * r3[3] * r3[3];
//if (iterations > 50)
// it = 0; // debug
if (disct > 1e-8)
if (disct > 0)
{
// 3-roots, find the largest P...
it = 0;
@ -4112,7 +4111,7 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
// accept a (possible) whobble in the curve...
// error_msg("No convergence when calculating P in Peng-Robinson.", STOP);
}
if (V_m < v1)
if (V_m < v1 && it < 40)
P = R_TK / (v1 - b_sum) - a_aa_sum / (v1 * (v1 + 2 * b_sum) - b2);
}
}
@ -4162,11 +4161,11 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
continue;
}
phase_ptr->pr_p = phase_ptr->fraction_x * P;
if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
{
phase_ptr->pr_phi = 1;
continue;
}
//if (phase_ptr->t_c == 0.0 || phase_ptr->p_c == 0.0)
//{
// phase_ptr->pr_phi = 1;
// continue;
//}
rz = P * V_m / R_TK;
A = a_aa_sum * P / (R_TK * R_TK);
B = b_sum * P / R_TK;
@ -4181,6 +4180,8 @@ calc_PR(std::vector<struct phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
// for initial equilibrations, adapt log_k of the gas phase...
if (state < REACTION)
{
rho_0 = calc_rho_0(TK - 273.15, P);
calc_dielectrics(TK - 273.15, P);
phase_ptr->lk = calc_lk_phase(phase_ptr, TK, P);
}
phase_ptr->pr_in = true;
@ -4280,6 +4281,7 @@ adjust_setup_pure_phases(void)
if (phase_ptr->p_c > 0 && phase_ptr->t_c > 0)
{
p = exp(si_org * LOG_10);
patm_x = p;
t = use.Get_solution_ptr()->Get_tc() + 273.15;
if (!phase_ptr->pr_in || p != phase_ptr->pr_p || t != phase_ptr->pr_tk)
{
@ -4685,6 +4687,7 @@ adjust_setup_solution(void)
if (phase_ptr->p_c > 0 && phase_ptr->t_c > 0)
{
p = exp(x[i]->si * LOG_10);
patm_x = p;
t = use.Get_solution_ptr()->Get_tc() + 273.15;
if (!phase_ptr->pr_in || p != phase_ptr->pr_p || t != phase_ptr->pr_tk)
{
@ -5747,7 +5750,8 @@ calc_vm(LDBLE tc, LDBLE pa)
s_x[i]->rxn_x->logk[vm_tc] += (s_x[i]->logk[vm3] + (s_x[i]->logk[vm4] + s_x[i]->logk[vm5] * tc) * tc) * mu_x;
/* for calculating delta_v of the reaction... */
s_search(s_x[i]->name)->logk[vm_tc] = s_x[i]->rxn_x->logk[vm_tc];
//s_search(s_x[i]->name)->logk[vm_tc] = s_x[i]->rxn_x->logk[vm_tc];
s_x[i]->logk[vm_tc] = s_x[i]->rxn_x->logk[vm_tc];
}
return OK;
}
@ -5761,11 +5765,11 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
* Calculates log k's for all species and pure_phases
*/
#ifndef PHREEQC2
if (same_model == TRUE && same_temperature == TRUE && same_pressure == TRUE && same_mu)
return (OK);
if (tc == current_tc && pa == current_pa && ((fabs(mu_x - current_mu) < 1e-3 * mu_x) || !mu_terms_in_logk))
return OK;
#else
if (same_model == TRUE && same_temperature == TRUE && same_pressure == TRUE)
return (OK);
if (tc == current_tc && pa == current_pa)
return OK;
#endif
int i;
LDBLE tempk = tc + 273.15;
@ -5780,14 +5784,15 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
#endif
calc_vm(tc, pa);
mu_terms_in_logk = false;
for (i = 0; i < count_s_x; i++)
{
if (s_x[i]->rxn_x->logk[vm_tc])
/* calculate delta_v for the reaction... */
s_x[i]->rxn_x->logk[delta_v] = calc_delta_v(s_x[i]->rxn_x, false);
if (same_model && same_temperature && s_x[i]->rxn_x->logk[delta_v] == 0)
if (tc == current_tc && s_x[i]->rxn_x->logk[delta_v] == 0)
continue;
mu_terms_in_logk = true;
s_x[i]->lk = k_calc(s_x[i]->rxn_x->logk, tempk, pa * PASCAL_PER_ATM);
}
/*
@ -5805,6 +5810,8 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
phases[i]->rxn_x->logk[delta_v] = calc_delta_v(phases[i]->rxn_x, true) -
(phases[i]->logk[vm0] + (phases[i]->logk[vm1] + phases[i]->logk[vm2] * tc) * tc -
phases[i]->logk[kappa] * (pa - 1));
if (phases[i]->rxn_x->logk[delta_v])
mu_terms_in_logk = true;
phases[i]->lk = k_calc(phases[i]->rxn_x->logk, tempk, pa * PASCAL_PER_ATM);
#endif
}
@ -5823,6 +5830,11 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
}
}
}
current_tc = tc;
current_pa = pa;
current_mu = mu_x;
return (OK);
}
@ -6036,6 +6048,12 @@ save_model(void)
last_model.count_surface_charge = 0;
last_model.surface_charge = NULL;
}
current_tc = NAN;
current_pa = NAN;
current_mu = NAN;
mu_terms_in_logk = true;
return (OK);
}
@ -6051,32 +6069,8 @@ check_same_model(void)
if (last_model.force_prep == TRUE)
{
last_model.force_prep = FALSE;
same_temperature = FALSE;
same_pressure = FALSE;
return (FALSE);
}
/*
* Check temperature
*/
if (fabs(tc_x - last_model.temperature) < 0.05)
{
same_temperature = TRUE;
}
else
{
same_temperature = FALSE;
}
/*
* Check pressure
*/
if (fabs(patm_x - last_model.pressure) < 0.01)
{
same_pressure = TRUE;
}
else
{
same_pressure = FALSE;
}
/*
* Check master species
*/

358
phreeqc/read.cpp
View File

@ -4997,364 +4997,6 @@ read_solution(void)
return (return_value);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_solution(void)
/* ---------------------------------------------------------------------- */
{
/*
* Reads solution data
*
* Arguments:
* none
*
* Returns:
* KEYWORD if keyword encountered, input_error may be incremented if
* a keyword is encountered in an unexpected position
* EOF if eof encountered while reading mass balance concentrations
* ERROR if error occurred reading data
*
*/
int i, j, n, l;
int n_user, n_user_end;
int default_pe, alk;
int count_isotopes;
int max_mass_balance, count_mass_balance;
char *ptr, *ptr1;
char *description;
char token[MAX_LENGTH], token1[MAX_LENGTH];
int return_value, opt;
char *next_char;
const char *opt_list[] = {
"temp", /* 0 */
"temperature", /* 1 */
"dens", /* 2 */
"density", /* 3 */
"units", /* 4 */
"redox", /* 5 */
"ph", /* 6 */
"pe", /* 7 */
"unit", /* 8 */
"isotope", /* 9 */
"water", /* 10 */
"press", /* 11 */
"pressure" /* 12 */
};
int count_opt_list = 13;
/*
* Read solution number and description
*/
ptr = line;
read_number_description(ptr, &n_user, &n_user_end, &description);
/*
* Malloc space for solution data
*/
if (solution_bsearch(n_user, &n, FALSE) != NULL)
{
solution_free(solution[n]);
}
else
{
n = count_solution++;
if (count_solution >= max_solution)
{
space((void **) ((void *) &(solution)), count_solution,
&max_solution, sizeof(struct solution *));
}
}
solution[n] = solution_alloc();
solution[n]->n_user = n_user;
solution[n]->n_user_end = n_user_end;
if (use.Get_solution_in() == FALSE)
{
use.Set_solution_in(true);
use.Set_n_solution_user(n_user);
}
max_mass_balance = MAX_MASS_BALANCE;
/*
* Set default ph, temp, density, pe, units
*/
solution[n]->description = description;
solution[n]->tc = 25.0;
solution[n]->patm = 1;
solution[n]->ph = 7.0;
solution[n]->density = 1.0;
solution[n]->solution_pe = 4.0;
solution[n]->mass_water = 1.0;
solution[n]->ah2o = 1.0;
solution[n]->mu = 1e-7;
solution[n]->cb = 0.0;
solution[n]->count_master_activity = 0;
default_pe = 0;
solution[n]->units = string_hsave("mMol/kgw");
solution[n]->totals[0].description = NULL;
count_mass_balance = 0;
count_isotopes = 0;
/*
* Read concentration data
*/
return_value = UNKNOWN;
for (;;)
{
opt = get_option(opt_list, count_opt_list, &next_char);
if (opt == OPTION_DEFAULT)
{
ptr = next_char;
if (copy_token(token, &ptr, &l) == DIGIT)
{
opt = 9;
}
}
switch (opt)
{
case OPTION_EOF: /* end of file */
return_value = EOF;
break;
case OPTION_KEYWORD: /* keyword */
return_value = KEYWORD;
break;
case OPTION_ERROR:
input_error++;
error_msg("Unknown input in SOLUTION keyword.", CONTINUE);
error_msg(line_save, CONTINUE);
break;
case 0: /* temperature */
case 1:
if (sscanf(next_char, SCANFORMAT, &(solution[n]->tc)) != 1)
{
solution[n]->tc = 25;
}
break;
case 2: /* density */
case 3:
sscanf(next_char, SCANFORMAT, &(solution[n]->density));
break;
case 4: /* units */
case 8: /* unit */
if (copy_token(token, &next_char, &l) == EMPTY)
break;
{
if (check_units(token, FALSE, FALSE, solution[n]->units, TRUE) ==
OK)
{
solution[n]->units = string_hsave(token);
}
else
{
input_error++;
}
}
break;
case 5: /* redox */
if (copy_token(token, &next_char, &l) == EMPTY)
break;
if (parse_couple(token) == OK)
{
default_pe = pe_data_store(&(solution[n]->pe), token);
}
else
{
input_error++;
}
break;
case 6: /* ph */
if (read_conc(n, count_mass_balance, line) == ERROR)
{
input_error++;
break;
}
solution[n]->ph =
solution[n]->totals[count_mass_balance].input_conc;
if (solution[n]->totals[count_mass_balance].equation_name == NULL)
{
break;
}
solution[n]->totals[count_mass_balance].description =
string_hsave("H(1)");
count_mass_balance++;
break;
case 7: /* pe */
if (read_conc(n, count_mass_balance, line) == ERROR)
{
input_error++;
break;
}
solution[n]->solution_pe =
solution[n]->totals[count_mass_balance].input_conc;
if (solution[n]->totals[count_mass_balance].equation_name == NULL)
{
break;
}
solution[n]->totals[count_mass_balance].description =
string_hsave("E");
count_mass_balance++;
break;
case 9: /* isotope */
if (copy_token(token, &next_char, &l) != DIGIT)
{
input_error++;
error_string = sformatf( "Expected isotope name to"
" begin with an isotopic number.");
error_msg(error_string, CONTINUE);
continue;
}
solution[n]->isotopes =
(struct isotope *) PHRQ_realloc(solution[n]->isotopes,
(size_t) (count_isotopes +
1) *
sizeof(struct isotope));
if (solution[n]->isotopes == NULL)
malloc_error();
/* read and save element name */
ptr1 = token;
get_num(&ptr1,
&(solution[n]->isotopes[count_isotopes].isotope_number));
if (ptr1[0] == '\0' || isupper((int) ptr1[0]) == FALSE)
{
error_msg("Expecting element name.", CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
return (ERROR);
}
solution[n]->isotopes[count_isotopes].elt_name =
string_hsave(ptr1);
/* read and store isotope ratio */
if (copy_token(token, &next_char, &l) != DIGIT)
{
input_error++;
error_string = sformatf(
"Expected numeric value for isotope ratio.");
error_msg(error_string, CONTINUE);
continue;
}
sscanf(token, SCANFORMAT,
&(solution[n]->isotopes[count_isotopes].ratio));
solution[n]->isotopes[count_isotopes].ratio_uncertainty = NAN;
/* read and store isotope ratio uncertainty */
if ((j = copy_token(token, &next_char, &l)) != EMPTY)
{
if (j != DIGIT)
{
input_error++;
error_string = sformatf(
"Expected numeric value for uncertainty in isotope ratio.");
error_msg(error_string, CONTINUE);
continue;
}
sscanf(token, SCANFORMAT,
&(solution[n]->isotopes[count_isotopes].
ratio_uncertainty));
}
count_isotopes++;
break;
case 10: /* water */
j = copy_token(token, &next_char, &l);
if (j == EMPTY)
{
solution[n]->mass_water = 1.0;
}
else if (j != DIGIT)
{
input_error++;
error_string = sformatf(
"Expected numeric value for mass of water in solution.");
error_msg(error_string, CONTINUE);
}
else
{
sscanf(token, SCANFORMAT, &dummy);
solution[n]->mass_water = (LDBLE) dummy;
}
break;
case 11: /* pressure */
case 12:
j = copy_token(token, &next_char, &l);
if (j == EMPTY)
{
solution[n]->patm = 1.0;
}
else if (j != DIGIT)
{
input_error++;
error_string = sformatf(
"Expected numeric value for pressure (atm) of solution.");
error_msg(error_string, CONTINUE);
}
else
{
sscanf(token, SCANFORMAT, &dummy);
solution[n]->patm = (LDBLE) dummy;
}
break;
case OPTION_DEFAULT:
/*
* Read concentration
*/
if (read_conc(n, count_mass_balance, line) == ERROR)
{
input_error++;
break;
}
count_mass_balance++;
break;
}
if (count_mass_balance + 1 >= max_mass_balance)
{
space((void **) ((void *) &(solution[n]->totals)),
count_mass_balance + 1, &max_mass_balance,
sizeof(struct conc));
}
if (return_value == EOF || return_value == KEYWORD)
break;
}
/*
* fix up default units and default pe
*/
for (i = 0; i < count_mass_balance; i++)
{
strcpy(token, solution[n]->totals[i].description);
str_tolower(token);
if (solution[n]->totals[i].units == NULL)
{
solution[n]->totals[i].units = solution[n]->units;
}
else
{
alk = FALSE;
if (strstr(token, "alk") == token)
alk = TRUE;
strcpy(token1, solution[n]->totals[i].units);
if (check_units(token1, alk, TRUE, solution[n]->units, TRUE) ==
ERROR)
{
input_error++;
}
else
{
solution[n]->totals[i].units = string_hsave(token1);
}
}
if (solution[n]->totals[i].n_pe < 0)
{
solution[n]->totals[i].n_pe = default_pe;
}
}
solution[n]->default_pe = default_pe;
/*
* Mark end of solution
*/
solution[n]->totals[count_mass_balance].description = NULL;
solution[n]->count_isotopes = count_isotopes;
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_species(void)

14
phreeqc/sit.cpp
View File

@ -504,7 +504,7 @@ set_sit(int initial)
#ifdef PHREEQC2
patm_x = solution_ptr->Get_patm();
#else
//patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
#endif
/*
* H+, e-, H2O
@ -763,6 +763,7 @@ jacobian_sit(void)
int i, j;
Restart:
int pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x);
if (full_pitzer == TRUE)
{
molalities(TRUE);
@ -800,6 +801,8 @@ Restart:
d2 = d1;
break;
case PITZER_GAMMA:
if (!full_pitzer)
continue;
x[i]->s->lg += d;
d2 = d;
break;
@ -833,8 +836,10 @@ Restart:
x[i]->moles += d2;
break;
case MU:
//continue;
d2 = d * mu_x;
mu_x += d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
break;
case PP:
@ -888,6 +893,7 @@ Restart:
break;
case MU:
mu_x -= d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
break;
case GAS_MOLES:
@ -1157,11 +1163,6 @@ check_gammas_sit(void)
{
converge = FALSE;
}
if (fabs(old_mu - mu_x) > 1e-3 * mu_x)
{
same_model = false;
k_temp(tc_x, patm_x);
}
t = pow((LDBLE) 10.0, s_h2o->la);
if ((pow((LDBLE) 10.0, s_h2o->la) - AW) > tol)
{
@ -1181,6 +1182,7 @@ gammas_sit()
int i, j;
LDBLE coef;
/* Initialize */
k_temp(tc_x, patm_x);
/*
* Calculate activity coefficients
*/

1
phreeqc/step.cpp
View File

@ -1015,7 +1015,6 @@ add_gas_phase(cxxGasPhase *gas_phase_ptr)
}
if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE && fabs(gas_phase_ptr->Get_total_p() - patm_x) > 0.01)
{
same_pressure = FALSE;
patm_x = gas_phase_ptr->Get_total_p();
k_temp(tc_x, patm_x);
}

32
phreeqc/transport.cpp
View File

@ -1439,8 +1439,8 @@ fill_spec(int l_cell_no)
continue;
}
dum2 = s_ptr->moles * dum; /* equivalent fraction */
sol_D[l_cell_no].spec[count_spec].name =
string_hsave(s_ptr->name);
sol_D[l_cell_no].spec[count_spec].name = s_ptr->name;
//string_hsave(s_ptr->name);
sol_D[l_cell_no].spec[count_spec].type = EX;
sol_D[l_cell_no].spec[count_spec].c = dum2;
sol_D[l_cell_no].spec[count_spec].lg = s_ptr->lg - log10(dum);
@ -1466,8 +1466,8 @@ fill_spec(int l_cell_no)
break;
}
/* copy its name and Dw and charge... */
sol_D[l_cell_no].spec[count_spec].aq_name =
string_hsave(s_ptr2->name);
sol_D[l_cell_no].spec[count_spec].aq_name = s_ptr2->name;
//string_hsave(s_ptr2->name);
sol_D[l_cell_no].spec[count_spec].z = s_ptr2->z;
if (s_ptr2->dw == 0)
sol_D[l_cell_no].spec[count_spec].Dwt =
@ -1484,7 +1484,8 @@ fill_spec(int l_cell_no)
lm = s_ptr->lm;
if (lm > MIN_LM)
{
sol_D[l_cell_no].spec[count_spec].name = string_hsave(s_ptr->name);
//sol_D[l_cell_no].spec[count_spec].name = string_hsave(s_ptr->name);
sol_D[l_cell_no].spec[count_spec].name = s_ptr->name;
sol_D[l_cell_no].spec[count_spec].type = AQ;
sol_D[l_cell_no].spec[count_spec].c =
s_ptr->moles / mass_water_aq_x;
@ -1847,7 +1848,8 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec)
}
if (l == count_m_s)
{
m_s[l].name = string_hsave(elt_list[k].elt->name);
//m_s[l].name = string_hsave(elt_list[k].elt->name);
m_s[l].name = elt_list[k].elt->name;
m_s[l].tot1 = elt_list[k].coef * l_J_ij[j].tot1;
m_s[l].tot2 = elt_list[k].coef * l_J_ij[j].tot2;
count_m_s++;
@ -2297,7 +2299,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
/* species 'name' is only in icell */
if (il_calcs && sol_D[icell].spec[i].type == EX)
{
J_ij_il[k_il].name = string_hsave(sol_D[icell].spec[i].name);
//J_ij_il[k_il].name = string_hsave(sol_D[icell].spec[i].name);
J_ij_il[k_il].name = sol_D[icell].spec[i].name;
V_M_il[k_il].D = sol_D[icell].spec[i].Dwt;
V_M_il[k_il].z = sol_D[icell].spec[i].z;
V_M_il[k_il].Dz = V_M_il[k_il].D * V_M_il[k_il].z;
@ -2312,7 +2315,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
else
{
J_ij[k].name = string_hsave(sol_D[icell].spec[i].name);
//J_ij[k].name = string_hsave(sol_D[icell].spec[i].name);
J_ij[k].name = sol_D[icell].spec[i].name;
V_M[k].D = sol_D[icell].spec[i].Dwt;
V_M[k].z = sol_D[icell].spec[i].z;
V_M[k].Dz = V_M[k].D * V_M[k].z;
@ -2376,7 +2380,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
/* species 'name' is only in jcell */
if (il_calcs && sol_D[jcell].spec[j].type == EX)
{
J_ij_il[k_il].name = string_hsave(sol_D[jcell].spec[j].name);
//J_ij_il[k_il].name = string_hsave(sol_D[jcell].spec[j].name);
J_ij_il[k_il].name = sol_D[jcell].spec[j].name;
V_M_il[k_il].D = sol_D[jcell].spec[j].Dwt;
V_M_il[k_il].z = sol_D[jcell].spec[j].z;
V_M_il[k_il].Dz = V_M_il[k_il].D * V_M_il[k_il].z;
@ -2391,7 +2396,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
else
{
J_ij[k].name = string_hsave(sol_D[jcell].spec[j].name);
//J_ij[k].name = string_hsave(sol_D[jcell].spec[j].name);
J_ij[k].name = sol_D[jcell].spec[j].name;
V_M[k].D = sol_D[jcell].spec[j].Dwt;
V_M[k].z = sol_D[jcell].spec[j].z;
V_M[k].Dz = V_M[k].D * V_M[k].z;
@ -2455,7 +2461,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
/* species 'name' is in both cells */
if (il_calcs && sol_D[icell].spec[i].type == EX)
{
J_ij_il[k_il].name = string_hsave(sol_D[icell].spec[i].name);
//J_ij_il[k_il].name = string_hsave(sol_D[icell].spec[i].name);
J_ij_il[k_il].name = sol_D[icell].spec[i].name;
if (sol_D[icell].spec[i].Dwt == 0
|| sol_D[jcell].spec[j].Dwt == 0)
V_M_il[k_il].D = 0.0;
@ -2478,7 +2485,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
else
{
J_ij[k].name = string_hsave(sol_D[icell].spec[i].name);
//J_ij[k].name = string_hsave(sol_D[icell].spec[i].name);
J_ij[k].name = sol_D[icell].spec[i].name;
if (sol_D[icell].spec[i].Dwt == 0
|| sol_D[jcell].spec[j].Dwt == 0)
V_M[k].D = 0.0;

12
phreeqc/utilities.cpp
View File

@ -121,14 +121,20 @@ calc_rho_0(LDBLE tc, LDBLE pa)
LDBLE p2 = -5.50294e-6 + tc * ( 1.07699e-7 + tc * (-2.05409e-9 + tc * ( 1.30238e-11 + tc * -3.20982E-14)));
/* The minimal pressure equals the saturation pressure... */
p_sat = exp(11.6702 - 3816.44 / (T - 46.13)) * ah2o_x;
ah2o_x0 = ah2o_x; // for updating rho in model(): compare with new ah2o_x
if (ah2o_x <= 1.0)
p_sat = exp(11.6702 - 3816.44 / (T - 46.13)) * ah2o_x;
else
p_sat = exp(11.6702 - 3816.44 / (T - 46.13));
//ah2o_x0 = ah2o_x; // for updating rho in model(): compare with new ah2o_x
if (pa < p_sat || (use.Get_solution_ptr() && use.Get_solution_ptr()->Get_patm() < p_sat))
{
pa = p_sat;
}
patm_x = pa;
if (!use.Get_gas_phase_in())
patm_x = pa;
rho_0 = p0 + pa * (p1 + pa * p2);
if (rho_0 < 0.01)
rho_0 = 0.01;
/* compressibility, d(ln(rho)) / d(P), 1/atm... */
kappa_0 = (p1 + 2 * pa * p2) / rho_0;