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

Runs most kinetics, kinsurf is a problem.

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/branches/ErrorHandling@6084 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2012-01-20 00:56:16 +00:00
parent 18906cdd2f
commit dd06c2eaf2
23 changed files with 3973 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
KineticsComp.cxx
View File

@ -27,12 +27,13 @@ PHRQ_base(io)
//
{
tol = 1e-8;
m = 0.0;
m0 = 0.0;
m = -1;
m0 = -1;
moles = 0.0;
initial_moles = 0;
namecoef.type = cxxNameDouble::ND_NAME_COEF;
}
#ifdef SKIP
cxxKineticsComp::cxxKineticsComp(struct kinetics_comp *kinetics_comp_ptr, PHRQ_io *io)
//
// constructor for cxxKineticsComp from struct kinetics_comp
@ -51,7 +52,7 @@ namecoef(kinetics_comp_ptr->list, kinetics_comp_ptr->count_list)
this->d_params.push_back(kinetics_comp_ptr->d_params[i]);
}
}
#endif
cxxKineticsComp::~cxxKineticsComp()
{
}
@ -118,8 +119,8 @@ cxxKineticsComp::dump_raw(std::ostream & s_oss, unsigned int indent) const
// Kinetics_Comp element and attributes
s_oss << indent0 << "-rate_name " << this->
rate_name << "\n";
//s_oss << indent0 << "-rate_name " << this->
// rate_name << "\n";
s_oss << indent1 << "-tol " << this->tol << "\n";
s_oss << indent1 << "-m " << this->m << "\n";
s_oss << indent1 << "-m0 " << this->m0 << "\n";
@ -160,7 +161,7 @@ cxxKineticsComp::read_raw(CParser & parser, bool check)
if (vopts.empty())
{
vopts.reserve(10);
vopts.push_back("rate_name"); // 0
vopts.push_back("rate_name_not_used"); // 0
vopts.push_back("tol"); // 1
vopts.push_back("m"); // 2
vopts.push_back("m0"); // 3
@ -205,7 +206,7 @@ cxxKineticsComp::read_raw(CParser & parser, bool check)
//parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE);
break;
case 0: // rate_name
case 0: // rate_name not used
if (!(parser.get_iss() >> str))
{
this->rate_name.clear();

24
KineticsComp.h
View File

@ -35,12 +35,20 @@ public:
this->rate_name.clear();
}
const cxxNameDouble &Get_namecoef(void) const {return namecoef;};
LDBLE Get_tol(void) const {return tol;};
LDBLE Get_m(void) const {return m;};
LDBLE Get_m0(void) const {return m0;};
LDBLE Get_moles(void) const {return moles;};
const std::vector < LDBLE > &Get_d_params(void) const {return d_params;};
cxxNameDouble &Get_namecoef(void) {return namecoef;}
void Set_namecoef(const cxxNameDouble nd) {namecoef = nd;}
LDBLE Get_tol(void) const {return tol;}
void Set_tol(LDBLE t) {tol = t;}
LDBLE Get_m(void) const {return m;}
void Set_m(LDBLE t) {m = t;}
LDBLE Get_m0(void) const {return m0;}
void Set_m0(LDBLE t) {m0 = t;}
LDBLE Get_moles(void) const {return moles;}
void Set_moles(LDBLE t) {moles = t;}
LDBLE Get_initial_moles(void) const {return initial_moles;}
void Set_initial_moles(LDBLE t) {initial_moles = t;}
std::vector < LDBLE > &Get_d_params(void) {return d_params;};
std::vector < std::string > &Get_c_params(void) {return c_params;};
#ifdef USE_MPI
void mpi_unpack(int *ints, int *ii, LDBLE *doubles, int *dd);
@ -56,7 +64,9 @@ public:
LDBLE m;
LDBLE m0;
LDBLE moles;
std::vector < LDBLE >d_params;
LDBLE initial_moles;
std::vector< LDBLE > d_params;
std::vector< std::string > c_params;
public:

1
PHRQ_io.cpp
View File

@ -75,6 +75,7 @@ output_msg(const char * str)
{
(*output_ostream) << str;
}
output_flush();
}
// ---------------------------------------------------------------------- */
// log ostream methods

22
Phreeqc.cpp
View File

@ -248,13 +248,23 @@ size_t Phreeqc::list_components(std::list<std::string> &list_c)
}
}
// kinetics
{
std::map<int, cxxKinetics>::iterator it = Rxn_kinetics_map.begin();
for (; it != Rxn_kinetics_map.end(); it++)
{
calc_dummy_kinetic_reaction_tally(&(it->second));
cxxKinetics entity = it->second;
accumulator.add_extensive(entity.Get_totals(), 1.0);
}
}
#ifdef SKIP
for (i = 0; i < count_kinetics; i++)
{
calc_dummy_kinetic_reaction_tally(&kinetics[i]);
cxxKinetics entity(&kinetics[i], phrq_io);
accumulator.add_extensive(entity.Get_totals(), 1.0);
}
#endif
// Put in all primaries
cxxNameDouble::iterator it;
for (it = accumulator.begin(); it != accumulator.end(); it++)
@ -328,7 +338,7 @@ void Phreeqc::init(void)
//max_exchange = MAX_PP_ASSEMBLAGE;
max_surface = MAX_PP_ASSEMBLAGE;
//max_gas_phase = MAX_PP_ASSEMBLAGE;
max_kinetics = MAX_PP_ASSEMBLAGE;
//max_kinetics = MAX_PP_ASSEMBLAGE;
//max_ss_assemblage = MAX_PP_ASSEMBLAGE;
max_elements = MAX_ELEMENTS;
@ -348,7 +358,7 @@ void Phreeqc::init(void)
//count_exchange = 0;
count_surface = 0;
//count_gas_phase = 0;
count_kinetics = 0;
//count_kinetics = 0;
//count_ss_assemblage = 0;
count_elements = 0;
@ -400,7 +410,7 @@ void Phreeqc::init(void)
//exchange = 0;
surface = 0;
//gas_phase = 0;
kinetics = 0;
//kinetics = 0;
//ss_assemblage = 0;
cell_data = 0;
elements = 0;
@ -554,7 +564,7 @@ void Phreeqc::init(void)
initial_total_time = 0;
rate_m = 0;
rate_m0 = 0;
rate_p = NULL;
//rate_p = NULL;
rate_time = 0;
rate_sim_time_start = 0;
rate_sim_time_end = 0;
@ -718,7 +728,7 @@ void Phreeqc::init(void)
//dbg_exchange = exchange;
dbg_surface = surface;
//dbg_pp_assemblage = pp_assemblage;
dbg_kinetics = kinetics;
//dbg_kinetics = kinetics;
//dbg_irrev = irrev;
//dbg_mix = mix;
dbg_master = master;

22
Phreeqc.h
View File

@ -390,8 +390,8 @@ public:
realtype uround, void *jac_data, long int *nfePtr,
N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
int calc_final_kinetic_reaction(struct kinetics *kinetics_ptr);
int calc_kinetic_reaction(struct kinetics *kinetics_ptr,
int calc_final_kinetic_reaction(cxxKinetics *kinetics_ptr);
int calc_kinetic_reaction(cxxKinetics *kinetics_ptr,
LDBLE time_step);
int rk_kinetics(int i, LDBLE kin_time, int use_mix, int nsaver,
LDBLE step_fraction);
@ -825,7 +825,7 @@ public:
//int add_exchange(struct exchange *exchange_ptr);
int add_exchange(cxxExchange *exchange_ptr);
int add_gas_phase(cxxGasPhase *gas_phase_ptr);
int add_kinetics(struct kinetics *kinetics_ptr);
int add_kinetics(cxxKinetics *kinetics_ptr);
//int add_mix(struct mix *mix_ptr);
int add_mix(cxxMix * mix_ptr);
//int add_pp_assemblage(struct pp_assemblage *pp_assemblage_ptr);
@ -872,6 +872,7 @@ public:
struct inverse *inverse_search(int n_user, int *n);
int inverse_sort(void);
public:
#ifdef SKIP
struct kinetics *kinetics_alloc(void);
struct kinetics *kinetics_bsearch(int k, int *n);
protected:
@ -895,6 +896,8 @@ protected:
int n_user_new);
struct kinetics *kinetics_search(int n_user, int *n, int print);
int kinetics_sort(void);
#endif
protected:
struct logk *logk_alloc(void);
int logk_copy2orig(struct logk *logk_ptr);
struct logk *logk_store(char *name, int replace_if_found);
@ -1050,8 +1053,8 @@ public:
// convert class to struct (structures.cpp)
//struct mix * cxxMix2mix(const cxxMix *mx);
struct kinetics *cxxKinetics2kinetics(const cxxKinetics * kin);
struct kinetics_comp * cxxKineticsComp2kinetics_comp(const std::list < cxxKineticsComp > * el);
//struct kinetics *cxxKinetics2kinetics(const cxxKinetics * kin);
//struct kinetics_comp * cxxKineticsComp2kinetics_comp(const std::list < cxxKineticsComp > * el);
//struct exchange * cxxExchange2exchange(const cxxExchange * ex);
//struct exch_comp * cxxExchComp2exch_comp(const std::map < std::string, cxxExchComp > * el);
//struct master * Get_exch_master(const cxxExchComp * ec);
@ -1084,7 +1087,7 @@ public:
/* tally.cpp */
void add_all_components_tally(void);
int build_tally_table(void);
int calc_dummy_kinetic_reaction_tally(struct kinetics *kinetics_ptr);
int calc_dummy_kinetic_reaction_tally(cxxKinetics *kinetics_ptr);
int diff_tally_table(void);
int extend_tally_table(void);
int free_tally_table(void);
@ -1293,11 +1296,13 @@ protected:
/* ----------------------------------------------------------------------
* Kinetics
* ---------------------------------------------------------------------- */
std::map<int, cxxKinetics> Rxn_kinetics_map;
#ifdef SKIP
struct kinetics *kinetics;
struct kinetics *dbg_kinetics;
int count_kinetics;
int max_kinetics;
#endif
int count_save_values;
struct save_values *save_values;
@ -1622,8 +1627,9 @@ protected:
struct rate *rates;
int count_rates;
LDBLE rate_m, rate_m0, *rate_p, rate_time, rate_sim_time_start,
LDBLE rate_m, rate_m0, /* *rate_p, */ rate_time, rate_sim_time_start,
rate_sim_time_end, rate_sim_time, rate_moles, initial_total_time;
std::vector<LDBLE> rate_p;
int count_rate_p;
/* ----------------------------------------------------------------------
* USER PRINT COMMANDS

2
Reaction.cxx
View File

@ -307,7 +307,7 @@ cxxReaction::read_raw(CParser & parser, const bool check)
}
}
int cxxReaction::
Get_actualSteps(void) const
Get_reaction_steps(void) const
{
if (equalIncrements)
{

2
Reaction.h
View File

@ -27,7 +27,7 @@ class cxxReaction:public cxxNumKeyword
cxxNameDouble &Get_reactantList(void) {return this->reactantList;}
std::vector < LDBLE > &Get_steps(void) {return this->steps;}
void Set_steps(std::vector<LDBLE> &v) {steps = v;}
int Get_actualSteps(void) const;
int Get_reaction_steps(void) const;
int Get_countSteps(void) const {return this->countSteps;}
void Set_countSteps(int i) {countSteps = i;}
bool Get_equalIncrements(void) const {return this->equalIncrements;}

64
ReadClass.cxx
View File

@ -207,6 +207,7 @@ read_surface_raw(void)
if (return_value == KEYWORD) echo_msg(sformatf( "\t%s\n", line));
return (return_value);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_kinetics_raw(void)
@ -298,6 +299,7 @@ read_kinetics_raw(void)
if (return_value == KEYWORD) echo_msg(sformatf( "\t%s\n", line));
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_dump(void)
@ -630,6 +632,7 @@ read_surface_modify(void)
if (return_value == OPTION_KEYWORD) echo_msg(sformatf( "\t%s\n", line));
return (return_value);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_kinetics_modify(void)
@ -718,6 +721,7 @@ read_kinetics_modify(void)
if (return_value == OPTION_KEYWORD) echo_msg(sformatf( "\t%s\n", line));
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
streamify_to_next_keyword(std::istringstream & lines)
@ -932,7 +936,23 @@ delete_entities(void)
}
}
}
// kineticss
// kinetics
if (delete_info.Get_kinetics().Get_defined())
{
if (delete_info.Get_kinetics().Get_numbers().size() == 0)
{
Rxn_kinetics_map.clear();
}
else
{
std::set < int >::iterator it;
for (it = delete_info.Get_kinetics().Get_numbers().begin(); it != delete_info.Get_kinetics().Get_numbers().end(); it++)
{
Rxn_kinetics_map.erase(*it);
}
}
}
#ifdef SKIP
if (delete_info.Get_kinetics().Get_defined())
{
if (delete_info.Get_kinetics().Get_numbers().size() == 0)
@ -954,6 +974,7 @@ delete_entities(void)
}
}
}
#endif
// mixes
if (delete_info.Get_mix().Get_defined())
{
@ -1029,9 +1050,9 @@ run_as_cells(void)
{
struct save save_data;
LDBLE kin_time;
int count_steps, use_mix, m;
int count_steps, use_mix;
char token[2 * MAX_LENGTH];
struct kinetics *kinetics_ptr;
//struct kinetics *kinetics_ptr;
state = REACTION;
if (run_info.Get_cells().Get_numbers().size() == 0 ||
@ -1069,14 +1090,14 @@ run_as_cells(void)
count_steps = 1;
if (use.Get_reaction_in() == TRUE && use.Get_reaction_ptr() != NULL)
{
int count = ((cxxReaction *) use.Get_reaction_ptr())->Get_actualSteps();
int count = ((cxxReaction *) use.Get_reaction_ptr())->Get_reaction_steps();
if (count > count_steps)
count_steps = count;
}
if (use.Get_kinetics_in() == TRUE && use.Get_kinetics_ptr() != NULL)
{
if (abs(use.Get_kinetics_ptr()->count_steps) > count_steps)
count_steps = abs(use.Get_kinetics_ptr()->count_steps);
if (use.Get_kinetics_ptr()->Get_reaction_steps() > count_steps)
count_steps = use.Get_kinetics_ptr()->Get_reaction_steps();
}
if (use.Get_temperature_in() == TRUE && use.Get_temperature_ptr() != NULL)
{
@ -1138,6 +1159,10 @@ run_as_cells(void)
// runner kin_time not defined
else
{
cxxKinetics *kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
kin_time = kinetics_ptr->Current_step((incremental_reactions==TRUE), reaction_step);
#ifdef SKIP
kinetics_ptr = kinetics_bsearch(-2, &m);
if (incremental_reactions == FALSE)
{
@ -1191,6 +1216,7 @@ run_as_cells(void)
}
}
}
#endif
}
}
if (incremental_reactions == FALSE ||
@ -1232,7 +1258,8 @@ run_as_cells(void)
memcpy(&save, &save_data, sizeof(struct save));
if (use.Get_kinetics_in() == TRUE)
{
kinetics_duplicate(-2, use.Get_n_kinetics_user());
//kinetics_duplicate(-2, use.Get_n_kinetics_user());
Utilities::Rxn_copy(Rxn_kinetics_map, -2, use.Get_n_kinetics_user());
}
saver();
}
@ -1387,6 +1414,27 @@ dump_ostream(std::ostream& os)
}
// kinetics
if (dump_info.Get_bool_kinetics())
{
if (dump_info.Get_kinetics().size() == 0)
{
Utilities::Rxn_dump_raw(Rxn_kinetics_map, os, 0);
}
else
{
std::set < int >::iterator it;
for (it = dump_info.Get_kinetics().begin(); it != dump_info.Get_kinetics().end(); it++)
{
cxxKinetics *p = Utilities::Rxn_find(Rxn_kinetics_map, *it);
if (p != NULL)
{
p->dump_raw(os, 0);
}
}
}
}
#ifdef SKIP
if (dump_info.Get_bool_kinetics())
{
if (dump_info.Get_kinetics().size() == 0)
@ -1411,7 +1459,7 @@ dump_ostream(std::ostream& os)
}
}
}
#endif
// mix
if (dump_info.Get_bool_mix())
{

7
Use.h
View File

@ -8,6 +8,7 @@ class cxxGasPhase;
class cxxPressure;
class cxxTemperature;
class cxxSSassemblage;
class cxxKinetics;
class cxxUse
{
@ -79,7 +80,7 @@ public:
cxxMix * Get_mix_ptr(void) const {return this->mix_ptr;}
cxxReaction * Get_reaction_ptr(void) const {return this->reaction_ptr;}
cxxExchange * Get_exchange_ptr(void) const {return this->exchange_ptr;}
struct kinetics * Get_kinetics_ptr(void) const {return this->kinetics_ptr;}
cxxKinetics * Get_kinetics_ptr(void) const {return this->kinetics_ptr;}
struct surface * Get_surface_ptr(void) const {return this->surface_ptr;}
cxxPressure * Get_pressure_ptr(void) const {return this->pressure_ptr;}
cxxTemperature * Get_temperature_ptr(void) const {return this->temperature_ptr;}
@ -92,7 +93,7 @@ public:
void Set_mix_ptr(cxxMix * p) {this->mix_ptr = p;}
void Set_reaction_ptr(cxxReaction * p) {this->reaction_ptr = p;}
void Set_exchange_ptr(cxxExchange * p) {this->exchange_ptr = p;}
void Set_kinetics_ptr(struct kinetics * p) {this->kinetics_ptr = p;}
void Set_kinetics_ptr(cxxKinetics * p) {this->kinetics_ptr = p;}
void Set_surface_ptr(struct surface * p) {this->surface_ptr = p;}
void Set_pressure_ptr(cxxPressure * p) {this->pressure_ptr = p;}
void Set_temperature_ptr(cxxTemperature * p) {this->temperature_ptr = p;}
@ -124,7 +125,7 @@ protected:
bool kinetics_in;
int n_kinetics_user;
struct kinetics *kinetics_ptr;
cxxKinetics *kinetics_ptr;
bool surface_in;
int n_surface_user;

454
cxxKinetics.cxx
View File

@ -35,9 +35,10 @@ cxxKinetics::cxxKinetics(PHRQ_io *io)
cvode_steps = 100;
cvode_order = 5;
totals.type = cxxNameDouble::ND_ELT_MOLES;
equal_steps = 0;
equalIncrements = false;
count = 0;
}
#ifdef SKIP
cxxKinetics::cxxKinetics(struct kinetics *kinetics_ptr, PHRQ_io *io)
//
// constructor for cxxKinetics from struct kinetics
@ -80,6 +81,7 @@ totals(kinetics_ptr->totals)
this->steps.push_back(kinetics_ptr->steps[i]);
}
}
#endif
cxxKinetics::cxxKinetics(const std::map < int, cxxKinetics > &entities,
cxxMix & mix, int l_n_user, PHRQ_io *io):
cxxNumKeyword(io)
@ -186,26 +188,32 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out) co
s_oss << "-cvode_order " << this->cvode_order << "\n";
// kineticsComps structures
for (std::list < cxxKineticsComp >::const_iterator it =
kineticsComps.begin(); it != kineticsComps.end(); ++it)
//for (std::list < cxxKineticsComp >::const_iterator it =
// kineticsComps.begin(); it != kineticsComps.end(); ++it)
for (size_t k = 0; k < this->kinetics_comps.size(); k++)
{
s_oss << indent1;
s_oss << "-component" << "\n";
(*it).dump_raw(s_oss, indent + 2);
s_oss << "-component " << this->kinetics_comps[k].Get_rate_name() << "\n";
this->kinetics_comps[k].dump_raw(s_oss, indent + 2);
}
// totals
s_oss << indent1;
s_oss << "-totals " << "\n";
s_oss << "-totals " << "\n";
this->totals.dump_raw(s_oss, indent + 2);
// equal_steps
s_oss << indent1;
s_oss << "-equal_steps " << this->equal_steps << "\n";
s_oss << "-equalIncrements " << this->equalIncrements << "\n";
// equal_steps
s_oss << indent1;
s_oss << "-count " << this->count << "\n";
// steps
s_oss << indent1;
s_oss << "-steps " << "\n";
s_oss << "-steps " << "\n";
{
int i = 0;
s_oss << indent2;
@ -224,7 +232,295 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent, int * n_out) co
}
return;
}
void
cxxKinetics::read_raw(CParser & parser, bool check)
{
LDBLE d;
static std::vector < std::string > vopts;
if (vopts.empty())
{
vopts.reserve(15);
vopts.push_back("step_divide");
vopts.push_back("rk");
vopts.push_back("bad_step_max");
vopts.push_back("use_cvode");
vopts.push_back("component");
vopts.push_back("totals");
vopts.push_back("steps");
vopts.push_back("cvode_steps");
vopts.push_back("cvode_order");
vopts.push_back("equalIncrements");
vopts.push_back("count");
}
std::istream::pos_type ptr;
std::istream::pos_type next_char;
std::string token;
int opt_save;
bool useLastLine(false);
std::vector < LDBLE > temp_steps;
// Read kinetics number and description
this->read_number_description(parser);
opt_save = CParser::OPT_ERROR;
bool step_divide_defined(false);
bool rk_defined(false);
bool bad_step_max_defined(false);
bool use_cvode_defined(false);
bool cvode_steps_defined(false);
bool cvode_order_defined(false);
bool steps_defined(false);
bool equalIncrements_defined(false);
bool count_defined(false);
for (;;)
{
int opt;
if (useLastLine == false)
{
opt = parser.get_option(vopts, next_char);
}
else
{
opt = parser.getOptionFromLastLine(vopts, next_char, true);
}
if (opt == CParser::OPT_DEFAULT)
{
opt = opt_save;
}
useLastLine = false;
switch (opt)
{
case CParser::OPT_EOF:
break;
case CParser::OPT_KEYWORD:
break;
case CParser::OPT_DEFAULT:
case CParser::OPT_ERROR:
opt = CParser::OPT_EOF;
parser.error_msg("Unknown input in KINETICS_COMP_RAW keyword.",
PHRQ_io::OT_CONTINUE);
parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE);
break;
case 0: // step_divide
if (!(parser.get_iss() >> this->step_divide))
{
this->step_divide = 1.0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for step_divide.",
PHRQ_io::OT_CONTINUE);
}
step_divide_defined = true;
break;
case 1: // rk
if (!(parser.get_iss() >> this->rk))
{
this->rk = 3;
parser.incr_input_error();
parser.error_msg("Expected integer value for rk.",
PHRQ_io::OT_CONTINUE);
}
rk_defined = true;
break;
case 2: // bad_step_max
if (!(parser.get_iss() >> this->bad_step_max))
{
this->bad_step_max = 500;
parser.incr_input_error();
parser.error_msg("Expected integer value for bad_step_max.",
PHRQ_io::OT_CONTINUE);
}
bad_step_max_defined = true;
break;
case 3: // use_cvode
if (!(parser.get_iss() >> this->use_cvode))
{
this->use_cvode = false;
parser.incr_input_error();
parser.error_msg("Expected boolean value for use_cvode.",
PHRQ_io::OT_CONTINUE);
}
use_cvode_defined = true;
break;
case 4: // component
{
std::string str;
if (!(parser.get_iss() >> str))
{
parser.incr_input_error();
parser.error_msg("Expected string value for component name.",
PHRQ_io::OT_CONTINUE);
}
else
{
cxxKineticsComp temp_comp(this->io);
temp_comp.Set_rate_name(str.c_str());
cxxKineticsComp *comp_ptr = this->Find(str);
if (comp_ptr)
{
temp_comp = *comp_ptr;
}
temp_comp.read_raw(parser, false);
if (comp_ptr)
{
for (size_t j = 0; j < this->kinetics_comps.size(); j++)
{
if (Utilities::strcmp_nocase(this->kinetics_comps[j].Get_rate_name().c_str(), str.c_str()) == 0)
{
this->kinetics_comps[j] = temp_comp;
}
}
}
else
{
this->kinetics_comps.push_back(temp_comp);
}
useLastLine = true;
}
}
break;
case 5: // totals
if (this->totals.read_raw(parser, next_char) !=
CParser::PARSER_OK)
{
parser.incr_input_error();
parser.
error_msg
("Expected element name and molality for KineticsComp totals.",
PHRQ_io::OT_CONTINUE);
}
opt_save = 5;
break;
case 6: // steps
while (parser.copy_token(token, next_char) == CParser::TT_DIGIT)
{
//sscanf(token.c_str(), "%lf", &d);
//this->steps.push_back(d);
std::istringstream iss(token);
if (!(iss >> d))
{
parser.incr_input_error();
parser.error_msg("Expected numeric value for steps.",
PHRQ_io::OT_CONTINUE);
}
else
{
temp_steps.push_back(d);
steps_defined = true;
}
}
opt_save = 6;
break;
case 7: // cvode_steps
if (!(parser.get_iss() >> this->cvode_steps))
{
this->cvode_steps = 100;
parser.incr_input_error();
parser.error_msg("Expected integer value for cvode_steps.",
PHRQ_io::OT_CONTINUE);
}
cvode_steps_defined = true;
break;
case 8: // cvode_order
if (!(parser.get_iss() >> this->cvode_order))
{
this->cvode_order = 5;
parser.incr_input_error();
parser.error_msg("Expected integer value for cvode_order.",
PHRQ_io::OT_CONTINUE);
}
cvode_order_defined = true;
break;
case 9: // equalIncrements
if (!(parser.get_iss() >> this->equalIncrements))
{
this->use_cvode = false;
parser.incr_input_error();
parser.error_msg("Expected boolean value for equalIncrements.",
PHRQ_io::OT_CONTINUE);
}
equalIncrements_defined = true;
break;
case 10: // count
if (!(parser.get_iss() >> this->count))
{
this->count = 0;
parser.incr_input_error();
parser.error_msg("Expected integer value for count.",
PHRQ_io::OT_CONTINUE);
}
count_defined = true;
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
}
if (steps_defined)
{
this->steps = temp_steps;
}
if (check)
{
// members that must be defined
if (step_divide_defined == false)
{
parser.incr_input_error();
parser.error_msg("Step_divide not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
if (rk_defined == false)
{
parser.incr_input_error();
parser.error_msg("Rk not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
if (bad_step_max_defined == false)
{
parser.incr_input_error();
parser.error_msg("Bad_step_max not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
if (use_cvode_defined == false)
{
parser.incr_input_error();
parser.error_msg("Use_cvode not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
if (cvode_steps_defined == false)
{
parser.incr_input_error();
parser.error_msg("Cvode_steps not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
if (cvode_order_defined == false)
{
parser.incr_input_error();
parser.error_msg("Cvode_order not defined for KINETICS_RAW input.",
PHRQ_io::OT_CONTINUE);
}
// disable to remain backward compatible with PHAST restart files?
//if (equal_steps_defined == false)
//{
// parser.incr_input_error();
// parser.error_msg("Equal_steps not defined for KINETICS_RAW input.",
// PHRQ_io::OT_CONTINUE);
//}
}
}
#ifdef SKIP
void
cxxKinetics::read_raw(CParser & parser, bool check)
{
@ -357,11 +653,12 @@ cxxKinetics::read_raw(CParser & parser, bool check)
reread.set_echo_file(CParser::EO_NONE);
reread.set_echo_stream(CParser::EO_NONE);
std::list < cxxKineticsComp >::iterator kit;
//std::list < cxxKineticsComp >::iterator kit;
bool found = false;
for (kit = this->kineticsComps.begin(); kit != this->kineticsComps.end(); kit++)
size_t k;
for (k = 0; k < this->kineticsComps.size(); k++)
{
if (kit->Get_rate_name() == ec.Get_rate_name())
if (this->kineticsComps[k].Get_rate_name() == ec.Get_rate_name())
{
found = true;
break;
@ -369,7 +666,7 @@ cxxKinetics::read_raw(CParser & parser, bool check)
}
if (found)
{
kit->read_raw(reread, false);
this->kineticsComps[k].read_raw(reread, false);
}
else
{
@ -523,7 +820,7 @@ cxxKinetics::read_raw(CParser & parser, bool check)
//}
}
}
#endif
#ifdef USE_MPI
void
cxxKinetics::mpi_pack(std::vector < int >&ints,
@ -603,6 +900,50 @@ cxxKinetics::add(const cxxKinetics & addee, LDBLE extensive)
if (extensive == 0.0)
return;
//std::map < std::string, cxxKineticsComp> kineticsComps;
for (size_t i_add = 0; i_add < addee.kinetics_comps.size(); i_add++)
{
bool found(false);
size_t i;
for (i = 0; i < this->kinetics_comps.size(); i++)
{
if (this->kinetics_comps[i].Get_rate_name() == addee.kinetics_comps[i_add].Get_rate_name())
{
found = true;
break;
}
}
if (found)
{
this->kinetics_comps[i].add(addee.kinetics_comps[i_add], extensive);
}
else
{
cxxKineticsComp entity = addee.kinetics_comps[i_add];
entity.multiply(extensive);
this->kinetics_comps.push_back(entity);
}
}
this->steps = addee.steps;
this->step_divide = addee.step_divide;
this->rk = addee.rk;
this->bad_step_max = addee.bad_step_max;
this->use_cvode = addee.use_cvode;
this->cvode_steps = addee.cvode_steps;
this->cvode_order = addee.cvode_order;
this->equalIncrements = addee.equalIncrements;
this->count = addee.count;
}
#ifdef SKIP
void
cxxKinetics::add(const cxxKinetics & addee, LDBLE extensive)
//
// Add to existing ppassemblage to "this" ppassemblage
//
{
if (extensive == 0.0)
return;
//std::map < std::string, cxxKineticsComp> kineticsComps;
for (size_t i = 0; i < kineticsComps.size(); i++)
for (std::list < cxxKineticsComp >::const_iterator itadd =
addee.kineticsComps.begin(); itadd != addee.kineticsComps.end();
++itadd)
@ -645,7 +986,84 @@ cxxKinetics::add(const cxxKinetics & addee, LDBLE extensive)
this->cvode_order = addee.cvode_order;
this->equal_steps = addee.equal_steps;
}
//cxxNameDouble & cxxKinetics::get_totals(void)
//{
// return this->totals;
//}
#endif
cxxKineticsComp * cxxKinetics::
Find(const std::string &s)
{
for (size_t i = 0; i < this->kinetics_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->kinetics_comps[i].Get_rate_name().c_str(), s.c_str()) == 0)
{
return & (this->kinetics_comps[i]);
}
}
return NULL;
}
int cxxKinetics::
Get_reaction_steps(void) const
{
if(equalIncrements)
{
return count;
}
return ((int) steps.size());
}
LDBLE cxxKinetics::
Current_step(bool incremental_reactions, int reaction_step) const
{
if (this->steps.size() == 0)
return 1;
LDBLE kin_time = 1;
if (!incremental_reactions)
{
if (!this->equalIncrements)
{
if (reaction_step > (int) this->steps.size())
{
kin_time = this->steps[this->steps.size() - 1];
}
else
{
kin_time = this->steps[reaction_step - 1];
}
}
else
{
if (reaction_step > this->count)
{
kin_time = this->steps[0];
}
else
{
kin_time = (LDBLE) reaction_step * this->steps[0] / ((LDBLE) (this->count));
}
}
}
else
{
/* incremental reactions */
if (!this->equalIncrements)
{
if (reaction_step > (int) this->steps.size())
{
kin_time = this->steps[this->steps.size() - 1];
}
else
{
kin_time = this->steps[reaction_step - 1];
}
}
else
{
if (reaction_step > this->count)
{
kin_time = 0;
}
else
{
kin_time = this->steps[0] / ((LDBLE) (this->count));
}
}
}
return kin_time;
}

40
cxxKinetics.h
View File

@ -17,7 +17,7 @@ class cxxKinetics:public cxxNumKeyword
public:
cxxKinetics(PHRQ_io *io=NULL);
cxxKinetics(struct kinetics *, PHRQ_io *io=NULL);
//cxxKinetics(struct kinetics *, PHRQ_io *io=NULL);
cxxKinetics(const std::map < int, cxxKinetics > &entity_map, cxxMix & mx,
int n_user, PHRQ_io *io=NULL);
~cxxKinetics();
@ -31,17 +31,29 @@ class cxxKinetics:public cxxNumKeyword
//bool Get_related_phases(void) const;
//bool Get_related_rate(void) const;
const std::vector < LDBLE > &Get_steps(void) const {return steps;};
LDBLE Get_step_divide(void) const {return step_divide;};
std::vector < LDBLE > &Get_steps(void) {return steps;}
LDBLE Get_step_divide(void) const {return step_divide;}
void Set_step_divide(LDBLE t) {step_divide = t;}
int Get_rk(void) const {return rk;};
int Get_bad_step_max(void) const {return bad_step_max;};
bool Get_use_cvode(void) const {return use_cvode;};
int Get_cvode_steps(void) const {return cvode_steps;};
int Get_cvode_order(void) const {return cvode_order;};
const std::list < cxxKineticsComp > &Get_kineticsComps(void) const {return kineticsComps;};
const cxxNameDouble & Get_totals(void) const {return this->totals;};
int Get_equal_steps(void) const {return equal_steps;};
void Set_rk(int t) {rk = t;}
int Get_bad_step_max(void) const {return bad_step_max;}
void Set_bad_step_max(int t) {bad_step_max = t;}
bool Get_use_cvode(void) const {return use_cvode;}
void Set_use_cvode(bool tf) {use_cvode = tf;}
int Get_cvode_steps(void) const {return cvode_steps;}
void Set_cvode_steps(int t) {cvode_steps = t;}
int Get_cvode_order(void) const {return cvode_order;}
void Set_cvode_order(int t) {cvode_order = t;}
std::vector < cxxKineticsComp > &Get_kinetics_comps(void) {return kinetics_comps;}
cxxNameDouble & Get_totals(void) {return this->totals;}
void Set_totals(const cxxNameDouble &nd) {totals = nd;}
bool Get_equalIncrements(void) const {return equalIncrements;}
void Set_equalIncrements(bool tf) {equalIncrements = tf;}
int Get_count(void) const {return count;}
void Set_count(int i) {count = i;}
int Get_reaction_steps(void) const;
cxxKineticsComp * Find(const std::string &str);
LDBLE Current_step(const bool incremental_reactions, const int reaction_step) const;
#ifdef USE_MPI
void mpi_unpack(int *ints, int *ii, LDBLE *doubles, int *dd);
@ -51,9 +63,10 @@ class cxxKinetics:public cxxNumKeyword
void add(const cxxKinetics & addee, LDBLE extensive);
protected:
std::list < cxxKineticsComp > kineticsComps;
std::vector < cxxKineticsComp > kinetics_comps;
std::vector < LDBLE >steps;
int equal_steps;
int count;
bool equalIncrements;
cxxNameDouble totals;
LDBLE step_divide;
int rk;
@ -61,7 +74,6 @@ class cxxKinetics:public cxxNumKeyword
bool use_cvode;
int cvode_steps;
int cvode_order;
};
#endif // !defined(CXXKINETICS_H_INCLUDED)

5
phreeqc/advection.cpp
View File

@ -1,5 +1,7 @@
#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
@ -36,7 +38,8 @@ advection(void)
{
for (i = 1; i <= count_ad_cells; i++)
{
if (kinetics_bsearch(i, &n) != NULL)
//if (kinetics_bsearch(i, &n) != NULL)
if (Utilities::Rxn_find(Rxn_kinetics_map, i) == NULL)
{
input_error++;
error_string = sformatf(

26
phreeqc/basicsubs.cpp
View File

@ -7,6 +7,7 @@
#include "../GasPhase.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
activity(const char *species_name)
@ -974,7 +975,30 @@ get_calculate_value(const char *name)
}
return (calculate_value_ptr->value);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
kinetics_moles(const char *kinetics_name)
/* ---------------------------------------------------------------------- */
{
if (use.Get_kinetics_in() == FALSE || use.Get_kinetics_ptr() == NULL)
return (0);
for (size_t i = 0; i < use.Get_kinetics_ptr()->Get_kinetics_comps().size(); i++)
{
cxxKineticsComp *kinetics_comp_ptr = &(use.Get_kinetics_ptr()->Get_kinetics_comps()[i]);
if (strcmp_nocase
(kinetics_comp_ptr->Get_rate_name().c_str(), kinetics_name) == 0)
{
return (kinetics_comp_ptr->Get_m());
}
}
error_string = sformatf( "No data for rate %s in KINETICS keyword.",
kinetics_name);
warning_msg(error_string);
return (0);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
kinetics_moles(const char *kinetics_name)
@ -998,7 +1022,7 @@ kinetics_moles(const char *kinetics_name)
warning_msg(error_string);
return (0);
}
#endif
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
log_activity(const char *species_name)

3
phreeqc/global_structures.h
View File

@ -409,6 +409,7 @@ struct Charge_Group
LDBLE z;
LDBLE eq;
};
#ifdef SKIP
/* ----------------------------------------------------------------------
* Kinetics
* ---------------------------------------------------------------------- */
@ -446,7 +447,7 @@ struct kinetics_comp
int count_d_params;
LDBLE *d_params;
};
#endif
/*----------------------------------------------------------------------
* Save
*---------------------------------------------------------------------- */

1737
phreeqc/kinetics.cpp
View File

File diff suppressed because it is too large Load Diff

124
phreeqc/mainsubs.cpp
View File

@ -6,11 +6,12 @@
#include "PBasic.h"
#include "Temperature.h"
#include "Exchange.h"
#include "ExchComp.h"
//#include "ExchComp.h"
#include "GasPhase.h"
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
#if defined(WINDOWS) || defined(_WINDOWS)
#include <windows.h>
@ -71,7 +72,7 @@ initialize(void)
//max_exchange = MAX_PP_ASSEMBLAGE;
max_surface = MAX_PP_ASSEMBLAGE;
//max_gas_phase = MAX_PP_ASSEMBLAGE;
max_kinetics = MAX_PP_ASSEMBLAGE;
//max_kinetics = MAX_PP_ASSEMBLAGE;
//max_ss_assemblage = MAX_PP_ASSEMBLAGE;
max_elements = MAX_ELEMENTS;
@ -91,7 +92,7 @@ initialize(void)
//count_exchange = 0;
count_surface = 0;
//count_gas_phase = 0;
count_kinetics = 0;
//count_kinetics = 0;
//count_ss_assemblage = 0;
count_elements = 0;
@ -158,8 +159,8 @@ initialize(void)
//space((void **) ((void *) &gas_phase), INIT, &max_gas_phase,
// sizeof(struct gas_phase));
space((void **) ((void *) &kinetics), INIT, &max_kinetics,
sizeof(struct kinetics));
//space((void **) ((void *) &kinetics), INIT, &max_kinetics,
// sizeof(struct kinetics));
//space((void **) ((void *) &ss_assemblage), INIT, &max_ss_assemblage,
// sizeof(struct ss_assemblage));
@ -382,7 +383,7 @@ initialize(void)
initial_total_time = 0;
rate_m = 0;
rate_m0 = 0;
rate_p = NULL;
//rate_p = NULL;
rate_time = 0;
rate_sim_time_start = 0;
rate_sim_time_end = 0;
@ -594,7 +595,7 @@ initialize(void)
//dbg_exchange = exchange;
dbg_surface = surface;
//dbg_pp_assemblage = pp_assemblage;
dbg_kinetics = kinetics;
//dbg_kinetics = kinetics;
//dbg_irrev = irrev;
//dbg_mix = mix;
dbg_master = master;
@ -916,10 +917,7 @@ set_use(void)
*/
if (use.Get_kinetics_in() == TRUE)
{
int n_kinetics;
use.Set_kinetics_ptr(
kinetics_bsearch(use.Get_n_kinetics_user(), &n_kinetics));
//use.Set_n_kinetics(n_kinetics);
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user()));
if (use.Get_kinetics_ptr() == NULL)
{
error_string = sformatf( "Kinetics %d not found.",
@ -1400,11 +1398,11 @@ reactions(void)
* mixture,
* or irreversible reaction.
*/
int count_steps, use_mix, m;
int count_steps, use_mix;
char token[2 * MAX_LENGTH];
struct save save_data;
LDBLE kin_time;
struct kinetics *kinetics_ptr;
cxxKinetics *kinetics_ptr;
state = REACTION;
/* last_model.force_prep = TRUE; */
@ -1418,13 +1416,15 @@ reactions(void)
if (use.Get_reaction_in() == TRUE && use.Get_reaction_ptr() != NULL)
{
cxxReaction *reaction_ptr = (cxxReaction *) use.Get_reaction_ptr();
if (reaction_ptr->Get_actualSteps() > count_steps)
count_steps = reaction_ptr->Get_actualSteps();
if (reaction_ptr->Get_reaction_steps() > count_steps)
count_steps = reaction_ptr->Get_reaction_steps();
}
if (use.Get_kinetics_in() == TRUE && use.Get_kinetics_ptr() != NULL)
{
if (abs(use.Get_kinetics_ptr()->count_steps) > count_steps)
count_steps = abs(use.Get_kinetics_ptr()->count_steps);
//if (abs(use.Get_kinetics_ptr()->count_steps) > count_steps)
// count_steps = abs(use.Get_kinetics_ptr()->count_steps);
if (use.Get_kinetics_ptr()->Get_reaction_steps() > count_steps)
count_steps = use.Get_kinetics_ptr()->Get_reaction_steps();
}
if (use.Get_temperature_in() == TRUE && use.Get_temperature_ptr() != NULL)
{
@ -1466,9 +1466,13 @@ reactions(void)
* Determine time step for kinetics
*/
kin_time = 0.0;
if (use.Get_kinetics_in() == TRUE)
{
kinetics_ptr = kinetics_bsearch(-2, &m);
//kinetics_ptr = kinetics_bsearch(-2, &m);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
kin_time = kinetics_ptr->Current_step((incremental_reactions==TRUE), reaction_step);
#ifdef SKIP
if (incremental_reactions == FALSE)
{
if (kinetics_ptr->count_steps > 0)
@ -1528,6 +1532,7 @@ reactions(void)
}
}
}
#endif
}
if (incremental_reactions == FALSE ||
(incremental_reactions == TRUE && reaction_step == 1))
@ -1566,10 +1571,16 @@ reactions(void)
* save end of reaction
*/
memcpy(&save, &save_data, sizeof(struct save));
if (use.Get_kinetics_in() == TRUE)
{
Utilities::Rxn_copy(Rxn_kinetics_map, -2, use.Get_n_kinetics_user());
}
#ifdef SKIP
if (use.Get_kinetics_in() == TRUE)
{
kinetics_duplicate(-2, use.Get_n_kinetics_user());
}
#endif
saver();
/* free_model_allocs(); */
@ -1651,12 +1662,31 @@ saver(void)
n = save.n_ss_assemblage_user;
xss_assemblage_save(n);
Utilities::Rxn_copies(Rxn_ss_assemblage_map, save.n_ss_assemblage_user, save.n_ss_assemblage_user_end);
//for (i = save.n_ss_assemblage_user + 1;
// i <= save.n_ss_assemblage_user_end; i++)
//{
// ss_assemblage_duplicate(n, i);
//}
}
if (save.kinetics == TRUE && use.Get_kinetics_in() == TRUE
/*&& use.Get_kinetics_ptr() != NULL */)
{
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
//use.Set_kinetics_ptr(kinetics_bsearch(use.Get_n_kinetics_user(), &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user()));
}
else if (use.Get_kinetics_in() != FALSE)
{
//use.Set_kinetics_ptr(kinetics_bsearch(-2, &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, -2));
}
if (use.Get_kinetics_ptr() != NULL)
{
n = use.Get_kinetics_ptr()->Get_n_user();
for (i = save.n_kinetics_user; i <= save.n_kinetics_user_end; i++)
{
//kinetics_duplicate(n, i);
Utilities::Rxn_copy(Rxn_kinetics_map, n, i);
}
}
}
#ifdef SKIP
if (save.kinetics == TRUE && use.Get_kinetics_in() == TRUE
/*&& use.Get_kinetics_ptr() != NULL */)
{
@ -1677,6 +1707,7 @@ saver(void)
}
}
}
#endif
return (OK);
}
/* ---------------------------------------------------------------------- */
@ -2538,6 +2569,19 @@ copy_use(int i)
/*
* Find kinetics
*/
if (use.Get_kinetics_in() == TRUE)
{
//kinetics_duplicate(use.Get_n_kinetics_user(), i);
Utilities::Rxn_copy(Rxn_kinetics_map, use.Get_n_kinetics_user(), i);
save.kinetics = TRUE;
save.n_kinetics_user = i;
save.n_kinetics_user_end = i;
}
else
{
save.kinetics = FALSE;
}
#ifdef SKIP
if (use.Get_kinetics_in() == TRUE)
{
kinetics_duplicate(use.Get_n_kinetics_user(), i);
@ -2549,6 +2593,7 @@ copy_use(int i)
{
save.kinetics = FALSE;
}
#endif
/*
* Find surface
*/
@ -2683,7 +2728,7 @@ step_save_surf(int n_user)
*
* input: n_user is user solution number of target
*/
int i, j, k, n, m;
int i, j, k, n;
struct surface *surface_ptr;
/*
* Malloc space for solution data
@ -2719,7 +2764,29 @@ step_save_surf(int n_user)
/*
* Update grams
*/
/*if (surface_ptr->edl == TRUE && surface_ptr->related_rate == TRUE && use.Get_kinetics_ptr() != NULL) { */
if ((surface_ptr->type == DDL || surface_ptr->type == CD_MUSIC)
&& surface_ptr->related_rate == TRUE && use.Get_kinetics_ptr() != NULL)
{
for (j = 0; j < surface_ptr->count_comps; j++)
{
if (surface_ptr->comps[j].rate_name != NULL)
{
cxxKinetics *kinetics_ptr = use.Get_kinetics_ptr();
for (size_t m = 0; m < kinetics_ptr->Get_kinetics_comps().size(); m++)
{
cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[m]);
if (strcmp_nocase
(kinetics_comp_ptr->Get_rate_name().c_str(),
surface_ptr->comps[j].rate_name) != 0)
continue;
surface_ptr->charge[surface_ptr->comps[j].charge].grams =
kinetics_comp_ptr->Get_m();
break;
}
}
}
}
#ifdef SKIP
if ((surface_ptr->type == DDL || surface_ptr->type == CD_MUSIC)
&& surface_ptr->related_rate == TRUE && use.Get_kinetics_ptr() != NULL)
{
@ -2740,6 +2807,7 @@ step_save_surf(int n_user)
}
}
}
#endif
return (OK);
}
@ -2977,14 +3045,14 @@ copy_entities(void)
{
for (j = 0; j < copy_kinetics.count; j++)
{
if (kinetics_bsearch(copy_kinetics.n_user[j], &n) != NULL)
if (Utilities::Rxn_find(Rxn_kinetics_map, copy_kinetics.n_user[j]) != NULL)
{
for (i = copy_kinetics.start[j]; i <= copy_kinetics.end[j];
i++)
{
if (i == copy_kinetics.n_user[j])
continue;
kinetics_duplicate(copy_kinetics.n_user[j], i);
Utilities::Rxn_copy(Rxn_kinetics_map, copy_kinetics.n_user[j], i);
}
}
else
@ -3001,7 +3069,6 @@ copy_entities(void)
{
for (j = 0; j < copy_ss_assemblage.count; j++)
{
//if (ss_assemblage_bsearch(copy_ss_assemblage.n_user[j], &n) != NULL)
if (Utilities::Rxn_find(Rxn_ss_assemblage_map, copy_ss_assemblage.n_user[j]) != NULL)
{
for (i = copy_ss_assemblage.start[j];
@ -3009,7 +3076,6 @@ copy_entities(void)
{
if (i == copy_ss_assemblage.n_user[j])
continue;
//ss_assemblage_duplicate(copy_ss_assemblage.n_user[j], i);
Utilities::Rxn_copy(Rxn_ss_assemblage_map, copy_ss_assemblage.n_user[j], i);
}
}

403
phreeqc/print.cpp
View File

@ -10,7 +10,7 @@
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
array_print(LDBLE * array_l, int row_count, int column_count,
@ -94,7 +94,7 @@ int Phreeqc::
punch_all(void)
/* ---------------------------------------------------------------------- */
{
int i;
//int i;
//#ifndef PHREEQ98 /* if not PHREEQ98 use the standard declaration */
// if (pr.hdf == FALSE && (punch.in == FALSE || pr.punch == FALSE) && user_graph->commands == NULL)
@ -113,11 +113,13 @@ punch_all(void)
//#else /* if PHREEQ98 execute punch_user_graph first, that is, before punch.in and pr.punch is checked */
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
use.Set_kinetics_ptr(kinetics_bsearch(use.Get_n_kinetics_user(), &i));
//use.Set_kinetics_ptr(kinetics_bsearch(use.Get_n_kinetics_user(), &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user()));
}
else if (use.Get_kinetics_in() != FALSE)
{
use.Set_kinetics_ptr(kinetics_bsearch(-2, &i));
//use.Set_kinetics_ptr(kinetics_bsearch(-2, &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, -2));
}
#if defined PHREEQ98
if (pr.user_graph == TRUE)
@ -830,6 +832,218 @@ print_reaction(void)
print_kinetics(void)
/* ---------------------------------------------------------------------- */
{
/*
* prints kinetic reaction,
* should be called only on final kinetic step
*/
//int i, j;
LDBLE sim_time;
cxxKinetics *kinetics_ptr;
if (pr.kinetics == FALSE || pr.all == FALSE)
return (OK);
if (state < REACTION)
return (OK);
kinetics_ptr = NULL;
if (use.Get_kinetics_in() == TRUE)
{
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
//kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &i);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user());
}
else
{
//kinetics_ptr = kinetics_bsearch(-2, &i);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
}
}
if (kinetics_ptr == NULL)
return (OK);
/*
* determine time step
*/
if (state == TRANSPORT || state == PHAST)
{
kin_time_x = timest;
}
else if (state == ADVECTION)
{
kin_time_x = advection_kin_time;
}
sim_time = 0.;
if (run_info.Get_run_cells())
{
sim_time = rate_sim_time;
}
else
{
if (incremental_reactions == TRUE)
{
if (!kinetics_ptr->Get_equalIncrements())
{
for (int i = 0; i < reaction_step; i++)
{
if (i < (int) kinetics_ptr->Get_steps().size())
{
sim_time += kinetics_ptr->Get_steps()[i];
}
else
{
sim_time += kinetics_ptr->Get_steps().back();
}
}
}
else
{
if (reaction_step > kinetics_ptr->Get_count())
{
sim_time = kinetics_ptr->Get_steps().front();
}
else
{
sim_time =
reaction_step * kinetics_ptr->Get_steps().front() /
((LDBLE) (kinetics_ptr->Get_count()));
}
}
}
#ifdef SKIP
if (incremental_reactions == TRUE)
{
//if (kinetics_ptr->count_steps > 0)
if (!kinetics_ptr->Get_equalIncrements())
{
for (i = 0; i < reaction_step; i++)
{
if (i < kinetics_ptr->count_steps)
{
sim_time += kinetics_ptr->steps[i];
}
else
{
sim_time +=
kinetics_ptr->steps[kinetics_ptr->count_steps - 1];
}
}
}
else if (kinetics_ptr->count_steps < 0)
{
if (reaction_step > -kinetics_ptr->count_steps)
{
sim_time = kinetics_ptr->steps[0];
}
else
{
sim_time =
reaction_step * kinetics_ptr->steps[0] /
((LDBLE) (-kinetics_ptr->count_steps));
}
}
}
#endif
}
/*
* Print amount of reaction
*/
if (phast == FALSE)
{
output_msg(sformatf("Kinetics %d.\t%s\n\n",
use.Get_n_kinetics_user(), kinetics_ptr->Get_description().c_str()));
}
else
{
output_msg(sformatf("Kinetics.\n\n"));
}
/*
* Print reaction
*/
if (state == TRANSPORT)
{
output_msg(sformatf("\tTime: %g seconds\n",
(double) (initial_total_time + transport_step * timest)));
output_msg(sformatf("\tTime step: %g seconds\n\n",
(double) kin_time_x));
}
else if (state == ADVECTION)
{
output_msg(sformatf("\tTime: %g seconds\n",
(double) (initial_total_time +
advection_step * advection_kin_time)));
output_msg(sformatf("\tTime step: %g seconds\n\n",
(double) kin_time_x));
}
else if (state == PHAST)
{
output_msg(sformatf("\tTime: %g seconds\n",
(double) rate_sim_time_end));
output_msg(sformatf("\tTime step: %g seconds\n\n",
(double) kin_time_x));
}
else if (state == REACTION)
{
if (incremental_reactions == FALSE)
{
output_msg(sformatf("\tTime step: %g seconds\n\n",
(double) kin_time_x));
}
else
{
output_msg(sformatf(
"\tTime step: %g seconds (Incremented time: %g seconds)\n\n",
(double) kin_time_x, (double) sim_time));
}
}
output_msg(sformatf("\t%-15s%12s%12s %-15s%12s\n\n",
"Rate name", "Delta Moles", "Total Moles", "Reactant",
"Coefficient"));
//for (i = 0; i < kinetics_ptr->count_comps; i++)
for (size_t i = 0; i < kinetics_ptr->Get_kinetics_comps().size(); i++)
{
cxxKineticsComp *kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[i]);
if (state != TRANSPORT && state != PHAST)
{
output_msg(sformatf("\t%-15s%12.3e%12.3e",
kinetics_comp_ptr->Get_rate_name().c_str(),
(double) -kinetics_comp_ptr->Get_moles(),
(double) kinetics_comp_ptr->Get_m()));
}
else
{
output_msg(sformatf("\t%-15s%12.3e%12.3e",
kinetics_comp_ptr->Get_rate_name().c_str(),
(double) (kinetics_comp_ptr->Get_m() -
kinetics_comp_ptr->Get_initial_moles()),
(double) kinetics_comp_ptr->Get_m()));
}
//for (j = 0; j < kinetics_ptr->comps[i].count_list; j++)
cxxNameDouble::iterator it = kinetics_comp_ptr->Get_namecoef().begin();
for ( ; it != kinetics_comp_ptr->Get_namecoef().end(); it++)
{
std::string name = it->first;
LDBLE coef = it->second;
if (it == kinetics_comp_ptr->Get_namecoef().begin())
{
output_msg(sformatf(" %-15s%12g\n",
name.c_str(),
(double) coef));
}
else
{
output_msg(sformatf("\t%39s %-15s%12g\n", " ",
name.c_str(),
(double) coef));
}
}
}
output_msg(sformatf("\n"));
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
print_kinetics(void)
/* ---------------------------------------------------------------------- */
{
/*
* prints kinetic reaction,
* should be called only on final kinetic step
@ -995,7 +1209,7 @@ print_kinetics(void)
output_msg(sformatf("\n"));
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
print_master_reactions(void)
@ -2200,7 +2414,63 @@ print_totals(void)
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
print_user_print(void)
/* ---------------------------------------------------------------------- */
{
/*
* Print with user defined BASIC print routine
*/
cxxKinetics *kinetics_ptr;
char l_command[] = "run";
if (pr.user_print == FALSE || pr.all == FALSE)
return (OK);
if (user_print->commands == NULL)
return (OK);
kinetics_ptr = NULL;
if (use.Get_kinetics_in() == TRUE)
{
kinetics_ptr = use.Get_kinetics_ptr();
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
//use.Set_kinetics_ptr(kinetics_bsearch(use.Get_n_kinetics_user(), &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user()));
}
else
{
//use.Set_kinetics_ptr(kinetics_bsearch(-2, &i));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, -2));
}
}
print_centered("User print");
if (user_print->new_def == TRUE)
{
/* basic_renumber(user_print->commands, &user_print->linebase, &user_print->varbase, &user_print->loopbase); */
if (basic_compile
(user_print->commands, &user_print->linebase,
&user_print->varbase, &user_print->loopbase) != 0)
{
error_msg("Fatal Basic error in USER_PRINT.", STOP);
}
user_print->new_def = FALSE;
}
if (basic_run
(l_command, user_print->linebase, user_print->varbase,
user_print->loopbase) != 0)
{
error_msg("Fatal Basic error in USER_PRINT.", STOP);
}
output_msg(sformatf("\n"));
if (use.Get_kinetics_in() == TRUE)
{
use.Set_kinetics_ptr(kinetics_ptr);
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
print_user_print(void)
@ -2256,7 +2526,7 @@ print_user_print(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
print_using(void)
@ -2274,7 +2544,7 @@ print_using(void)
//struct ss_assemblage *ss_assemblage_ptr;
//struct gas_phase *gas_phase_ptr;
//struct irrev *irrev_ptr;
struct kinetics *kinetics_ptr;
//struct kinetics *kinetics_ptr;
int n;
if (pr.use == FALSE || pr.all == FALSE)
@ -2380,16 +2650,19 @@ print_using(void)
}
if (use.Get_kinetics_in() == TRUE)
{
cxxKinetics * kinetics_ptr;
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &n);
//kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &n);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user());
}
else
{
kinetics_ptr = kinetics_bsearch(-2, &n);
//kinetics_ptr = kinetics_bsearch(-2, &n);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
}
output_msg(sformatf("Using kinetics %d.\t%s\n",
use.Get_n_kinetics_user(), kinetics_ptr->description));
use.Get_n_kinetics_user(), kinetics_ptr->Get_description().c_str()));
}
output_msg(sformatf("\n"));
return (OK);
@ -2897,6 +3170,37 @@ punch_identifiers(void)
if (state == REACTION && incremental_reactions == TRUE
&& use.Get_kinetics_ptr() != NULL)
{
if (!use.Get_kinetics_ptr()->Get_equalIncrements())
{
reaction_time = 0.0;
for (i = 0; i < reaction_step; i++)
{
//if (i < use.Get_kinetics_ptr()->count_steps)
if (i < (int) use.Get_kinetics_ptr()->Get_steps().size())
{
reaction_time += use.Get_kinetics_ptr()->Get_steps()[i];
}
else
{
reaction_time +=
use.Get_kinetics_ptr()->Get_steps().back();
}
}
}
else
{
if (reaction_step > use.Get_kinetics_ptr()->Get_count())
{
reaction_time = use.Get_kinetics_ptr()->Get_steps().front();
}
else
{
reaction_time =
reaction_step * use.Get_kinetics_ptr()->Get_steps().front() /
((LDBLE) (use.Get_kinetics_ptr()->Get_count()));
}
}
#ifdef SKIP
if (use.Get_kinetics_ptr()->count_steps > 0)
{
reaction_time = 0.0;
@ -2926,6 +3230,7 @@ punch_identifiers(void)
((LDBLE) (-use.Get_kinetics_ptr()->count_steps));
}
}
#endif
}
if (state == REACTION)
{
@ -3140,7 +3445,80 @@ punch_saturation_indices(void)
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
punch_kinetics(void)
/* ---------------------------------------------------------------------- */
{
/*
* prints kinetic reaction,
* should be called only on final kinetic step
*/
//int i, j;
cxxKinetics *kinetics_ptr;
LDBLE moles, delta_moles;
kinetics_ptr = NULL;
if (use.Get_kinetics_in() == TRUE)
{
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
//kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &i);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user());
}
else
{
//kinetics_ptr = kinetics_bsearch(-2, &i);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
}
}
for (int i = 0; i < punch.count_kinetics; i++)
{
moles = 0.0;
delta_moles = 0.0;
if (kinetics_ptr != NULL)
{
for (size_t j = 0; j < kinetics_ptr->Get_kinetics_comps().size(); j++)
{
cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
if (strcmp_nocase
(punch.kinetics[i].name,
kinetics_comp_ptr->Get_rate_name().c_str()) == 0)
{
if (state != TRANSPORT && state != PHAST)
{
moles = kinetics_comp_ptr->Get_m();
delta_moles = - kinetics_comp_ptr->Get_moles();
}
else
{
moles = kinetics_comp_ptr->Get_m();
delta_moles =
kinetics_comp_ptr->Get_m() -
kinetics_comp_ptr->Get_initial_moles();
}
break;
}
}
}
if (punch.high_precision == FALSE)
{
fpunchf(sformatf("k_%s", punch.kinetics[i].name), "%12.4e\t",
(double) moles);
fpunchf(sformatf("dk_%s", punch.kinetics[i].name), "%12.4e\t",
(double) -delta_moles);
}
else
{
fpunchf(sformatf("k_%s", punch.kinetics[i].name), "%20.12e\t",
(double) moles);
fpunchf(sformatf("dk_%s", punch.kinetics[i].name), "%20.12e\t",
(double) -delta_moles);
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
punch_kinetics(void)
@ -3159,7 +3537,8 @@ punch_kinetics(void)
{
if (state == TRANSPORT || state == PHAST || state == ADVECTION)
{
kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &i);
//kinetics_ptr = kinetics_bsearch(use.Get_n_kinetics_user(), &i);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user());
}
else
{
@ -3211,7 +3590,7 @@ punch_kinetics(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
punch_user_punch(void)

497
phreeqc/read.cpp
View File

@ -14,8 +14,7 @@
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "SS.h"
#include "SScomp.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
@ -271,7 +270,8 @@ read_input(void)
Utilities::Rxn_read_raw(Rxn_pp_assemblage_map, this);
break;
case Keywords::KEY_KINETICS_RAW:
read_kinetics_raw();
//read_kinetics_raw();
Utilities::Rxn_read_raw(Rxn_kinetics_map, this);
break;
case Keywords::KEY_SOLID_SOLUTIONS_RAW:
Utilities::Rxn_read_raw(Rxn_ss_assemblage_map, this);
@ -310,7 +310,8 @@ read_input(void)
Utilities::Rxn_read_modify(Rxn_gas_phase_map, this);
break;
case Keywords::KEY_KINETICS_MODIFY:
read_kinetics_modify();
//read_kinetics_modify();
Utilities::Rxn_read_modify(Rxn_kinetics_map, this);
break;
case Keywords::KEY_DELETE:
read_delete();
@ -2115,7 +2116,490 @@ read_inv_phases(struct inverse *inverse_ptr, char *ptr)
inverse_ptr->count_phases++;
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_kinetics(void)
/* ---------------------------------------------------------------------- */
{
/*
* Reads kinetics 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
*
*/
/*
* Read kinetics
*/
char *ptr;
char *description;
std::string token;
int n_user, n_user_end;
LDBLE step;
int return_value, opt;
char *next_char;
const char *opt_list[] = {
"tol", /* 0 */
"m", /* 1 */
"m0", /* 2 */
"parms", /* 3 */
"formula", /* 4 */
"steps", /* 5 */
"step_divide", /* 6 */
"parameters", /* 7 */
"runge-kutta", /* 8 */
"runge_kutta", /* 9 */
"rk", /* 10 */
"bad_step_max", /* 11 */
"cvode", /* 12 */
"cvode_steps", /* 13 */
"cvode_order", /* 14 */
"time_steps" /* 15 */
};
int count_opt_list = 16;
/*
* Read kinetics number
*/
ptr = line;
read_number_description(ptr, &n_user, &n_user_end, &description);
cxxKinetics temp_kinetics;
temp_kinetics.Set_n_user(n_user);
temp_kinetics.Set_n_user_end(n_user_end);
temp_kinetics.Set_description(description);
description = (char *) free_check_null(description);
cxxKineticsComp *kinetics_comp_ptr = NULL;
/*
* Set use data to first read
*/
if (use.Get_kinetics_in() == FALSE)
{
use.Set_kinetics_in(true);
use.Set_n_kinetics_user(n_user);
}
/*
* Read kinetics data
*/
return_value = UNKNOWN;
for (;;)
{
opt = get_option(opt_list, count_opt_list, &next_char);
switch (opt)
{
case OPTION_EOF: /* end of file */
return_value = EOF;
break;
case OPTION_KEYWORD: /* keyword */
return_value = KEYWORD;
break;
case OPTION_DEFAULT: /* allocate space, read new name */
if (kinetics_comp_ptr)
{
temp_kinetics.Get_kinetics_comps().push_back(*kinetics_comp_ptr);
delete kinetics_comp_ptr;
}
kinetics_comp_ptr = new cxxKineticsComp;
ptr = line;
copy_token(token, &ptr);
kinetics_comp_ptr->Set_rate_name(token.c_str());
break;
case OPTION_ERROR:
input_error++;
error_msg("Unknown input in KINETICS keyword.", CONTINUE);
error_msg(line_save, CONTINUE);
break;
case 0: /* tolerance */
if (kinetics_comp_ptr == NULL)
{
error_string = sformatf( "No rate name has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
else
{
prev_next_char = next_char;
if (copy_token(token, &next_char) == DIGIT)
{
kinetics_comp_ptr->Set_tol(strtod(token.c_str(), &ptr));
}
else
{
error_string = sformatf(
"Expecting numerical value for tolerance, but found:\n %s",
prev_next_char);
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 1: /* m */
if (kinetics_comp_ptr == NULL)
{
error_string = sformatf( "No rate name has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
else
{
prev_next_char = next_char;
if (copy_token(token, &next_char) == DIGIT)
{
kinetics_comp_ptr->Set_m(strtod(token.c_str(), &ptr));
}
else
{
error_string = sformatf(
"Expecting numerical value for moles of reactant (m), but found:\n %s",
prev_next_char);
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 2: /* m0 */
if (kinetics_comp_ptr == NULL)
{
error_string = sformatf( "No rate name has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
else
{
prev_next_char = next_char;
if (copy_token(token, &next_char) == DIGIT)
{
kinetics_comp_ptr->Set_m0(strtod(token.c_str(), &ptr));
}
else
{
error_string = sformatf(
"Expecting numerical value for initial moles of reactant (m0), but found:\n %s",
prev_next_char);
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 3: /* parms */
case 7: /* parameters */
if (kinetics_comp_ptr == NULL)
{
error_string = sformatf( "No rate name has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
else
{
int j;
while ((j = copy_token(token, &next_char)) != EMPTY)
{
/*
* Store a LDBLE parameter
*/
if (j == DIGIT)
{
kinetics_comp_ptr->Get_d_params().push_back(strtod(token.c_str(), &ptr));
}
else
{
/*
* Store a character parameter
*/
kinetics_comp_ptr->Get_c_params().push_back(token);
}
}
}
break;
case 4: /* formula */
if (kinetics_comp_ptr == NULL)
{
error_string = sformatf( "No rate name has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
else
{
/*
* Store reactant name, default coefficient
*/
ptr = next_char;
bool have_name = false;
std::string name;
LDBLE coef;
while (copy_token(token, &ptr) != EMPTY)
{
coef = 1;
if (isalpha((int) token[0]) || (token[0] == '(')
|| (token[0] == '['))
{
if (have_name)
{
kinetics_comp_ptr->Get_namecoef().add(name.c_str(), coef);
}
name = token;
have_name = true;
}
else
{
if (!have_name)
{
error_string = sformatf( "No phase or chemical formula has been defined.");
error_msg(error_string, CONTINUE);
input_error++;
}
/*
* Store relative coefficient
*/
int j = sscanf(token.c_str(), SCANFORMAT, &coef);
if (j == 1)
{
kinetics_comp_ptr->Get_namecoef().add(name.c_str(), coef);
have_name = false;
}
else
{
error_msg
("Reading relative coefficient of reactant.",
CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
}
}
}
if (have_name)
{
kinetics_comp_ptr->Get_namecoef().add(name.c_str(), coef);
}
}
break;
case 5: /* steps */
case 15: /* time_steps */
/*
* Read one or more kinetics time increments
*/
{
int j;
while ((j = copy_token(token, &next_char)) == DIGIT)
{
/* Read next step increment(s) */
/* multiple, equal timesteps 15 aug. 2005 */
if (Utilities::replace("*", " ", token))
{
int k;
if (sscanf(token.c_str(), "%d" SCANFORMAT, &k, &step) == 2)
{
for (int i = 0; i < k; i++)
{
temp_kinetics.Get_steps().push_back(step);
}
}
else
{
input_error++;
error_msg
("Format error in multiple, equal KINETICS timesteps.\nCorrect is (for example): 20 4*10 2*5 3\n",
CONTINUE);
}
}
else
{
step = strtod(token.c_str(), &ptr);
temp_kinetics.Get_steps().push_back(step);
}
}
if (j == EMPTY)
break;
/*
* Read number of increments
*/
if (temp_kinetics.Get_steps().size() != 1)
{
error_msg
("To define equal time increments, only one total time should be defined.",
CONTINUE);
input_error++;
break;
}
temp_kinetics.Set_equalIncrements(true);
do
{
int i = 1;
j = sscanf(token.c_str(), "%d", &i);
if (j == 1)
{
temp_kinetics.Set_count(abs(i));
break;
}
else if (j == 1 && i < 0)
{
error_msg
("Expecting positive number for number of equal "
"time increments for kinetics.", CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
break;
}
}
while (copy_token(token, &next_char) != EMPTY);
}
break;
case 6: /* step_divide */
if (copy_token(token, &next_char) == DIGIT)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
temp_kinetics.Set_step_divide(dummy);
}
else
{
error_string = sformatf(
"Expecting numerical value for step_divide.");
error_msg(error_string, CONTINUE);
input_error++;
}
break;
case 8: /* runge-kutta */
case 9: /* runge_kutta */
case 10: /* rk */
{
int j = copy_token(token, &next_char);
if (j == DIGIT)
{
temp_kinetics.Set_rk((int) strtod(token.c_str(), &ptr));
}
else if (j == EMPTY)
{
}
else
{
error_string = sformatf(
"Expecting order for Runge-Kutta method.");
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 11: /* bad_step_max */
{
int j = copy_token(token, &next_char);
if (j == DIGIT)
{
temp_kinetics.Set_bad_step_max((int) strtod(token.c_str(), &ptr));
}
else if (j == EMPTY)
{
}
else
{
error_string = sformatf( "Expecting maximal bad steps number.");
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 12: /* cvode */
temp_kinetics.Set_use_cvode(get_true_false(next_char, TRUE) == TRUE);
break;
case 13: /* cvode_steps */
{
int j = copy_token(token, &next_char);
if (j == DIGIT)
{
temp_kinetics.Set_cvode_steps((int) strtod(token.c_str(), &ptr));
}
else if (j == EMPTY)
{
}
else
{
error_string = sformatf(
"Expecting maximum number of steps for one call to cvode.");
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
case 14: /* cvode_order */
{
int j = copy_token(token, &next_char);
if (j == DIGIT)
{
temp_kinetics.Set_cvode_order((int) strtod(token.c_str(), &ptr));
}
else if (j == EMPTY)
{
}
else
{
error_string = sformatf(
"Expecting number of terms (order) used in cvode (1-5).");
error_msg(error_string, CONTINUE);
input_error++;
}
}
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
}
// save last comp
if (kinetics_comp_ptr)
{
temp_kinetics.Get_kinetics_comps().push_back(*kinetics_comp_ptr);
delete kinetics_comp_ptr;
}
/*
* Default reactant
*/
for (size_t i = 0; i < temp_kinetics.Get_kinetics_comps().size(); i++)
{
cxxKineticsComp *kinetics_comp_ptr = &(temp_kinetics.Get_kinetics_comps()[i]);
if (kinetics_comp_ptr->Get_namecoef().size() == 0)
{
kinetics_comp_ptr->Get_namecoef().add(kinetics_comp_ptr->Get_rate_name().c_str(), 1.0);
}
}
/*
* Default 1 sec
*/
if (temp_kinetics.Get_steps().size() == 0)
{
temp_kinetics.Get_steps().push_back(1.0);
}
/*
* set defaults for moles
*/
for (size_t i = 0; i < temp_kinetics.Get_kinetics_comps().size(); i++)
{
cxxKineticsComp *kinetics_comp_ptr = &(temp_kinetics.Get_kinetics_comps()[i]);
if (kinetics_comp_ptr->Get_m0() < 0)
{
if (kinetics_comp_ptr->Get_m() < 0)
{
kinetics_comp_ptr->Set_m0(1);
}
else
{
kinetics_comp_ptr->Set_m0(kinetics_comp_ptr->Get_m());
}
}
if (kinetics_comp_ptr->Get_m() < 0)
{
kinetics_comp_ptr->Set_m(kinetics_comp_ptr->Get_m0());
}
}
Rxn_kinetics_map[n_user] = temp_kinetics;
return (return_value);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_kinetics(void)
@ -2659,7 +3143,7 @@ read_kinetics(void)
}
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
LDBLE * Phreeqc::
read_list_doubles(char **ptr, int *count_doubles)
@ -4354,7 +4838,7 @@ read_reaction_steps(cxxReaction *reaction_ptr)
/*
* Read number of equal increments, store as negative integer
*/
if (reaction_ptr->Get_actualSteps() != 1)
if (reaction_ptr->Get_reaction_steps() != 1)
{
error_msg
("To define equal increments, only one reaction increment should be defined.",
@ -9489,7 +9973,6 @@ read_solid_solutions(void)
copy_token(token, &ptr);
ss_ptr->Set_name(token);
ss_ptr->Set_total_moles(NAN);
std::cerr << "Read: " << ss_ptr->Get_total_moles() << "\n";
break;
}
if (return_value == EOF || return_value == KEYWORD)

62
phreeqc/step.cpp
View File

@ -15,6 +15,7 @@
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
@ -855,7 +856,7 @@ add_reaction(cxxReaction *reaction_ptr, int step_number, LDBLE step_fraction)
}
else if (reaction_ptr->Get_equalIncrements() && reaction_ptr->Get_steps().size()> 0)
{
if (step_number > (int) reaction_ptr->Get_actualSteps())
if (step_number > (int) reaction_ptr->Get_reaction_steps())
{
step_x = reaction_ptr->Get_steps()[0];
}
@ -863,7 +864,7 @@ add_reaction(cxxReaction *reaction_ptr, int step_number, LDBLE step_fraction)
{
step_x = reaction_ptr->Get_steps()[0] *
((LDBLE) step_number) /
((LDBLE) (reaction_ptr->Get_actualSteps()));
((LDBLE) (reaction_ptr->Get_reaction_steps()));
}
}
else
@ -876,9 +877,9 @@ add_reaction(cxxReaction *reaction_ptr, int step_number, LDBLE step_fraction)
/* Incremental reactions */
if (!reaction_ptr->Get_equalIncrements() && reaction_ptr->Get_steps().size()> 0)
{
if (step_number > (int) reaction_ptr->Get_actualSteps())
if (step_number > (int) reaction_ptr->Get_reaction_steps())
{
step_x = reaction_ptr->Get_steps()[reaction_ptr->Get_actualSteps() - 1];
step_x = reaction_ptr->Get_steps()[reaction_ptr->Get_reaction_steps() - 1];
}
else
{
@ -887,13 +888,13 @@ add_reaction(cxxReaction *reaction_ptr, int step_number, LDBLE step_fraction)
}
else if (reaction_ptr->Get_equalIncrements() && reaction_ptr->Get_steps().size()> 0)
{
if (step_number > (int) reaction_ptr->Get_actualSteps())
if (step_number > (int) reaction_ptr->Get_reaction_steps())
{
step_x = 0;
}
else
{
step_x = reaction_ptr->Get_steps()[0] / ((LDBLE) (reaction_ptr->Get_actualSteps()));
step_x = reaction_ptr->Get_steps()[0] / ((LDBLE) (reaction_ptr->Get_reaction_steps()));
}
}
else
@ -1188,6 +1189,54 @@ add_ss_assemblage(cxxSSassemblage *ss_assemblage_ptr)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
add_kinetics(cxxKinetics *kinetics_ptr)
/* ---------------------------------------------------------------------- */
{
/*
* Add kinetic reaction
*/
//int i;
struct master *master_ptr;
/*
* Add reaction to totals
*/
if (kinetics_ptr->Get_totals().size() == 0)
return (OK);
cxxNameDouble::iterator it = kinetics_ptr->Get_totals().begin();
for (; it != kinetics_ptr->Get_totals().end(); it++)
{
LDBLE coef = it->second;
struct element *elt_ptr = element_store(it->first.c_str());
if (elt_ptr == NULL || (master_ptr = elt_ptr->primary) == NULL)
{
input_error++;
error_string = sformatf(
"Element %s in kinetic reaction not found in database.",
it->first.c_str());
error_msg(error_string, STOP);
}
if (master_ptr->s == s_hplus)
{
total_h_x += coef;
}
else if (master_ptr->s == s_h2o)
{
total_o_x += coef;
}
else
{
master_ptr->total += coef;
if(master_ptr->total < 0.0)
{
fprintf(stderr, "Negative total\n");
}
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
add_kinetics(struct kinetics *kinetics_ptr)
/* ---------------------------------------------------------------------- */
{
@ -1227,6 +1276,7 @@ add_kinetics(struct kinetics *kinetics_ptr)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
gas_phase_check(cxxGasPhase *gas_phase_ptr)

76
phreeqc/structures.cpp
View File

@ -11,6 +11,7 @@
#include "PPassemblage.h"
#include "Use.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
@ -140,11 +141,14 @@ clean_up(void)
Rxn_gas_phase_map.clear();
/* kinetics */
Rxn_kinetics_map.clear();
#ifdef SKIP
for (j = 0; j < count_kinetics; j++)
{
kinetics_free(&kinetics[j]);
}
kinetics = (struct kinetics *) free_check_null(kinetics);
#endif
/* rates */
for (j = 0; j < count_rates; j++)
@ -362,7 +366,7 @@ clean_up(void)
//count_exchange = 0;
count_surface = 0;
//count_gas_phase = 0;
count_kinetics = 0;
//count_kinetics = 0;
//count_ss_assemblage = 0;
count_elements = 0;
@ -424,12 +428,14 @@ reinitialize(void)
Rxn_gas_phase_map.clear();
/* kinetics */
Rxn_kinetics_map.clear();
#ifdef SKIP
for (j = 0; j < count_kinetics; j++)
{
kinetics_free(&kinetics[j]);
}
count_kinetics = 0;
#endif
/* irreversible reactions */
Rxn_reaction_map.clear();
@ -970,6 +976,7 @@ inverse_sort(void)
}
return (OK);
}
#ifdef SKIP
/* **********************************************************************
*
* Routines related to structure "kinetics"
@ -1491,7 +1498,7 @@ kinetics_sort(void)
}
return (OK);
}
#endif
/* **********************************************************************
*
* Routines related to structure "master"
@ -4884,8 +4891,6 @@ system_duplicate(int i, int save_old)
if (solution_bsearch(i, &n, TRUE) != NULL)
solution_duplicate(i, save_old);
//if (pp_assemblage_bsearch(i, &n) != NULL)
// pp_assemblage_duplicate(i, save_old);
Utilities::Rxn_copy(Rxn_pp_assemblage_map, i, save_old);
Utilities::Rxn_copy(Rxn_exchange_map, i, save_old);
@ -4895,12 +4900,13 @@ system_duplicate(int i, int save_old)
Utilities::Rxn_copy(Rxn_gas_phase_map, i, save_old);
Utilities::Rxn_copy(Rxn_kinetics_map, i, save_old);
#ifdef SKIP
if (kinetics_bsearch(i, &n) != NULL)
kinetics_duplicate(i, save_old);
#endif
Utilities::Rxn_copy(Rxn_ss_assemblage_map, i, save_old);
//if (ss_assemblage_bsearch(i, &n) != NULL)
// ss_assemblage_duplicate(i, save_old);
return (OK);
}
@ -5174,7 +5180,7 @@ entity_exists(char *name, int n_user)
}
break;
case Kinetics: /* Kinetics */
if (kinetics_bsearch(n_user, &i) == NULL)
if (Utilities::Rxn_find(Rxn_kinetics_map, n_user) == NULL)
{
return_value = FALSE;
}
@ -5334,6 +5340,7 @@ copier_init(struct copier *copier_ptr)
(int *) PHRQ_malloc((size_t) (copier_ptr->max * sizeof(int)));
return (OK);
}
#ifdef SKIP
#include "../cxxKinetics.h"
struct kinetics * Phreeqc::
cxxKinetics2kinetics(const cxxKinetics * kin)
@ -5440,6 +5447,7 @@ cxxKineticsComp2kinetics_comp(const std::list < cxxKineticsComp > *el)
}
return (kinetics_comp_ptr);
}
#endif
#include "../Solution.h"
struct solution * Phreeqc::
cxxSolution2solution(const cxxSolution * sol)
@ -5928,6 +5936,15 @@ Use2cxxStorageBin(cxxStorageBin & sb)
sb.Set_SSassemblage(use.Get_n_ss_assemblage_user(), entity_ptr);
}
}
if (use.Get_kinetics_in())
{
cxxKinetics *entity_ptr = Utilities::Rxn_find(Rxn_kinetics_map, use.Get_n_kinetics_user());
if (entity_ptr != NULL)
{
sb.Set_Kinetics(use.Get_n_kinetics_user(), entity_ptr);
}
}
#ifdef SKIP
if (use.Get_kinetics_in())
{
struct kinetics *struct_entity = kinetics_bsearch(use.Get_n_kinetics_user(), &n);
@ -5937,6 +5954,7 @@ Use2cxxStorageBin(cxxStorageBin & sb)
sb.Set_Kinetics(use.Get_n_kinetics_user(), &entity);
}
}
#endif
if (use.Get_reaction_in())
{
cxxReaction *entity = Utilities::Rxn_find(Rxn_reaction_map, use.Get_n_reaction_user());
@ -5998,12 +6016,20 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb)
}
// Kinetics
{
std::map<int, cxxKinetics>::iterator it;
for (it = Rxn_kinetics_map.begin(); it != Rxn_kinetics_map.end(); it++)
{
sb.Set_Kinetics(it->second.Get_n_user(), &(it->second));
}
}
#ifdef SKIP
for (i = 0; i < count_kinetics; i++)
{
cxxKinetics entity(&kinetics[i], sb.Get_io());
sb.Set_Kinetics(kinetics[i].n_user, &entity );
}
#endif
// PPassemblages
{
std::map<int, cxxPPassemblage>::iterator it;
@ -6099,6 +6125,14 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb, int n)
}
// Kinetics
{
cxxKinetics *entity_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n);
if (entity_ptr != NULL)
{
sb.Set_Kinetics(n, entity_ptr);
}
}
#ifdef SKIP
{
if (kinetics_bsearch(n, &pos) != NULL)
{
@ -6107,7 +6141,7 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb, int n)
sb.Set_Kinetics(n, &ent );
}
}
#endif
// PPassemblages
{
cxxPPassemblage *entity_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n);
@ -6172,6 +6206,14 @@ cxxStorageBin2phreeqc(cxxStorageBin & sb, int n)
}
// Kinetics
{
std::map < int, cxxKinetics >::const_iterator it = sb.Get_Kinetics().find(n);
if (it != sb.Get_Kinetics().end())
{
Rxn_kinetics_map[n] = it->second;
}
}
#ifdef SKIP
{
std::map < int, cxxKinetics >::const_iterator it = sb.Get_Kinetics().find(n);
if (it != sb.Get_Kinetics().end())
@ -6182,7 +6224,7 @@ cxxStorageBin2phreeqc(cxxStorageBin & sb, int n)
kinetics_ptr = (struct kinetics *) free_check_null(kinetics_ptr);
}
}
#endif
// PPassemblages
{
std::map < int, cxxPPassemblage >::const_iterator it = sb.Get_PPassemblages().find(n);
@ -6284,6 +6326,14 @@ cxxStorageBin2phreeqc(cxxStorageBin & sb)
}
// Kinetics
{
std::map < int, cxxKinetics >::const_iterator it = sb.Get_Kinetics().begin();
for ( ; it != sb.Get_Kinetics().end(); it++)
{
Rxn_kinetics_map[it->first] = it->second;
}
}
#ifdef SKIP
{
std::map < int, cxxKinetics >::const_iterator it = sb.Get_Kinetics().begin();
for ( ; it != sb.Get_Kinetics().end(); it++)
@ -6294,7 +6344,7 @@ cxxStorageBin2phreeqc(cxxStorageBin & sb)
kinetics_ptr = (struct kinetics *) free_check_null(kinetics_ptr);
}
}
#endif
// PPassemblages
{
std::map < int, cxxPPassemblage >::const_iterator it = sb.Get_PPassemblages().begin();

217
phreeqc/tally.cpp
View File

@ -7,7 +7,7 @@
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/*
Calling sequence
@ -386,8 +386,8 @@ fill_tally_table(int *n_user, int index_conservative, int n_buffer)
struct surface *surface_ptr;
//struct ss_assemblage *ss_assemblage_ptr;
//struct gas_phase *gas_phase_ptr;
struct kinetics *kinetics_ptr;
struct kinetics_comp *kinetics_comp_ptr;
//struct kinetics *kinetics_ptr;
//struct kinetics_comp *kinetics_comp_ptr;
int found;
LDBLE moles;
char *ptr;
@ -614,31 +614,34 @@ fill_tally_table(int *n_user, int index_conservative, int n_buffer)
break;
}
case Kinetics:
/*
* fill in kinetics
*/
if (n_user[Kinetics] < 0)
break;
kinetics_ptr = kinetics_bsearch(n_user[Kinetics], &n);
if (kinetics_ptr == NULL)
break;
kinetics_comp_ptr = NULL;
for (j = 0; j < kinetics_ptr->count_comps; j++)
{
kinetics_comp_ptr = &kinetics_ptr->comps[j];
if (kinetics_comp_ptr->rate_name == tally_table[i].name)
/*
* fill in kinetics
*/
if (n_user[Kinetics] < 0)
break;
if (strcmp_nocase
(kinetics_comp_ptr->rate_name, tally_table[i].name) == 0)
//kinetics_ptr = kinetics_bsearch(n_user[Kinetics], &n);
cxxKinetics *kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n_user[Kinetics]);
if (kinetics_ptr == NULL)
break;
cxxKineticsComp * kinetics_comp_ptr = NULL;
for (j = 0; j < (int) kinetics_ptr->Get_kinetics_comps().size(); j++)
{
kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
if (string_hsave(kinetics_comp_ptr->Get_rate_name().c_str()) == tally_table[i].name)
break;
if (strcmp_nocase
(kinetics_comp_ptr->Get_rate_name().c_str(), tally_table[i].name) == 0)
break;
}
if (j >= (int) kinetics_ptr->Get_kinetics_comps().size())
break;
moles = kinetics_comp_ptr->Get_m();
count_elts = 0;
paren_count = 0;
add_elt_list(tally_table[i].formula, moles);
elt_list_to_tally_table(tally_table[i].total[n_buffer]);
}
if (j >= kinetics_ptr->count_comps)
break;
moles = kinetics_comp_ptr->m;
count_elts = 0;
paren_count = 0;
add_elt_list(tally_table[i].formula, moles);
elt_list_to_tally_table(tally_table[i].total[n_buffer]);
break;
case Mix:
break;
@ -707,7 +710,7 @@ build_tally_table(void)
* Also calculates a number greater than all user numbers and
* stores in global variable first_user_number.
*/
int i, j, k, l, n, p, save_print_use;
int j, k, l, n, p, save_print_use;
int count_tt_reaction, count_tt_exchange, count_tt_surface,
count_tt_gas_phase;
int count_tt_pure_phase, count_tt_ss_phase, count_tt_kinetics;
@ -716,8 +719,8 @@ build_tally_table(void)
//struct ss_assemblage *ss_assemblage_ptr;
//struct s_s *s_s_ptr;
//struct s_s_comp *s_s_comp_ptr;
struct kinetics *kinetics_ptr;
struct kinetics_comp *kinetics_comp_ptr;
//struct kinetics *kinetics_ptr;
//struct kinetics_comp *kinetics_comp_ptr;
struct phase *phase_ptr;
@ -936,6 +939,73 @@ build_tally_table(void)
* Add kinetic reactants
*/
count_tt_kinetics = 0;
if (Rxn_kinetics_map.size() > 0)
{
std::map<int, cxxKinetics>::iterator it;
for (it = Rxn_kinetics_map.begin(); it != Rxn_kinetics_map.end(); it++)
{
cxxKinetics *kinetics_ptr = &(it->second);
for (j = 0; j < (int) kinetics_ptr->Get_kinetics_comps().size(); j++)
{
cxxKineticsComp *kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
/*
* check if already in tally_table
*/
for (l = 1; l < count_tally_table_columns; l++)
{
if (tally_table[l].type == Kinetics &&
tally_table[l].name == string_hsave(kinetics_comp_ptr->Get_rate_name().c_str()))
break;
}
if (l < count_tally_table_columns)
continue;
/*
* Add to table
*/
count_tt_kinetics++;
n = count_tally_table_columns;
extend_tally_table();
tally_table[n].name = string_hsave(kinetics_comp_ptr->Get_rate_name().c_str());
tally_table[n].type = Kinetics;
/*
* get formula for kinetic component
*/
count_elts = 0;
paren_count = 0;
phase_ptr = NULL;
//if (kinetics_ptr->comps[j].count_list == 1)
if (kinetics_comp_ptr->Get_namecoef().size() == 1)
{
//strcpy(token, kinetics_ptr->comps[j].list[0].name);
strcpy(token, kinetics_comp_ptr->Get_namecoef().begin()->first.c_str());
phase_ptr = phase_bsearch(token, &p, FALSE);
}
if (phase_ptr != NULL)
{
add_elt_list(phase_ptr->next_elt, 1.0);
}
else
{
cxxNameDouble::iterator it = kinetics_comp_ptr->Get_namecoef().begin();
//for (k = 0; k < kinetics_ptr->comps[j].count_list; k++)
for ( ; it != kinetics_comp_ptr->Get_namecoef().end(); it++)
{
std::string name = it->first;
LDBLE coef = it->second;
char * temp_name = string_duplicate(name.c_str());
ptr = temp_name;
get_elts_in_species(&ptr, 1.0 * coef);
free_check_null(temp_name);
}
}
qsort(elt_list, (size_t) count_elts,
(size_t) sizeof(struct elt_list), elt_list_compare);
elt_list_combine();
tally_table[n].formula = elt_list_save();
}
}
}
#ifdef SKIP
if (count_kinetics > 0)
{
for (i = 0; i < count_kinetics; i++)
@ -998,6 +1068,7 @@ build_tally_table(void)
}
}
}
#endif
#ifdef SKIP
/*
* Debug print for table definition
@ -1117,21 +1188,88 @@ add_all_components_tally(void)
/*
* Add elements in kinetic reactions
*/
{
std::map<int, cxxKinetics>::iterator it = Rxn_kinetics_map.begin();
for ( ; it != Rxn_kinetics_map.end(); it++)
{
calc_dummy_kinetic_reaction_tally(&(it->second));
add_kinetics(&(it->second));
}
}
#ifdef SKIP
for (i = 0; i < count_kinetics; i++)
{
calc_dummy_kinetic_reaction_tally(&kinetics[i]);
add_kinetics(&kinetics[i]);
}
#endif
/*
* reset pr.use
*/
pr.use = save_print_use;
return;
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_dummy_kinetic_reaction_tally(struct kinetics *kinetics_ptr)
calc_dummy_kinetic_reaction_tally(cxxKinetics *kinetics_ptr)
/* ---------------------------------------------------------------------- */
{
/*
* Go through kinetic components and add positive amount of each reactant
*/
LDBLE coef;
//char token[MAX_LENGTH];
char *ptr;
struct phase *phase_ptr;
/*
* Go through list and generate list of elements and
* coefficient of elements in reaction
*/
//free_check_null(kinetics_ptr->totals);
kinetics_ptr->Get_totals().clear();
count_elts = 0;
paren_count = 0;
for (size_t i = 0; i < kinetics_ptr->Get_kinetics_comps().size(); i++)
{
cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[i]);
coef = 1.0;
/*
* Reactant is a pure phase, copy formula into token
*/
phase_ptr = NULL;
if (kinetics_comp_ptr->Get_namecoef().size() == 1)
{
std::string name = kinetics_comp_ptr->Get_namecoef().begin()->first;
int j;
phase_ptr = phase_bsearch(name.c_str(), &j, FALSE);
}
if (phase_ptr != NULL)
{
add_elt_list(phase_ptr->next_elt, coef);
}
else
{
cxxNameDouble::iterator it = kinetics_comp_ptr->Get_namecoef().begin();
for ( ; it != kinetics_comp_ptr->Get_namecoef().end(); it++)
{
std::string name = it->first;
char * temp_name = string_duplicate(name.c_str());
ptr = temp_name;
get_elts_in_species(&ptr, coef);
free_check_null(temp_name);
}
}
}
//kinetics_ptr->totals = elt_list_save();
kinetics_ptr->Set_totals(elt_list_NameDouble());
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_dummy_kinetic_reaction_tally(cxxKinetics *kinetics_ptr)
/* ---------------------------------------------------------------------- */
{
/*
@ -1146,7 +1284,8 @@ calc_dummy_kinetic_reaction_tally(struct kinetics *kinetics_ptr)
* Go through list and generate list of elements and
* coefficient of elements in reaction
*/
free_check_null(kinetics_ptr->totals);
//free_check_null(kinetics_ptr->totals);
kinetics_ptr->Get_totals().clear();
count_elts = 0;
paren_count = 0;
for (i = 0; i < kinetics_ptr->count_comps; i++)
@ -1181,7 +1320,7 @@ calc_dummy_kinetic_reaction_tally(struct kinetics *kinetics_ptr)
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
extend_tally_table(void)
@ -1256,6 +1395,21 @@ set_reaction_temperature(int n_user, LDBLE tc)
int Phreeqc::
set_kinetics_time(int n_user, LDBLE step)
/* ---------------------------------------------------------------------- */
{
cxxKinetics *kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n_user);
if (kinetics_ptr == NULL)
return (ERROR);
kinetics_ptr->Get_steps().clear();
kinetics_ptr->Get_steps().push_back(step);
kinetics_ptr->Set_equalIncrements(false);
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
set_kinetics_time(int n_user, LDBLE step)
/* ---------------------------------------------------------------------- */
{
int n;
struct kinetics *kinetics_ptr;
@ -1267,3 +1421,4 @@ set_kinetics_time(int n_user, LDBLE step)
kinetics_ptr->count_steps = 1;
return (OK);
}
#endif

370
phreeqc/tidy.cpp
View File

@ -5,7 +5,7 @@
#include "GasPhase.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
#define ZERO_TOL 1.0e-30
/* ---------------------------------------------------------------------- */
@ -13,7 +13,7 @@ int Phreeqc::
tidy_model(void)
/* ---------------------------------------------------------------------- */
{
int i, j;
int i;
int n_user, last;
int new_named_logk;
/*
@ -173,13 +173,13 @@ tidy_model(void)
phases[i]->pr_p = 0.0;
phases[i]->pr_phi = 1.0;
}
#ifdef SKIP
/* kinetics */
if (new_kinetics)
{
kinetics_sort();
}
#endif
/* named_log_k */
if (new_named_logk)
{
@ -289,18 +289,13 @@ tidy_model(void)
*/
if (new_kinetics)
{
for (i = 0; i < count_kinetics; i++)
std::map<int, cxxKinetics>::iterator it;
for (it = Rxn_kinetics_map.begin(); it != Rxn_kinetics_map.end(); it++)
{
if (kinetics[i].n_user_end > kinetics[i].n_user)
{
n_user = kinetics[i].n_user;
last = kinetics[i].n_user_end;
kinetics[i].n_user_end = kinetics[i].n_user;
for (j = n_user + 1; j <= last; j++)
{
kinetics_duplicate(n_user, j);
}
}
n_user = it->second.Get_n_user();
last = it->second.Get_n_user_end();
it->second.Set_n_user_end(n_user);
Utilities::Rxn_copies(Rxn_kinetics_map, n_user, last);
}
}
@ -2906,6 +2901,117 @@ tidy_isotopes(void)
int Phreeqc::
tidy_kin_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion
*/
{
//int k;
cxxKinetics *kinetics_ptr;
char *ptr;
LDBLE conc;
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
for ( ; it != Rxn_exchange_map.end(); it++)
{
cxxExchange * exchange_ptr = &(it->second);
if (!exchange_ptr->Get_new_def())
continue;
if (exchange_ptr->Get_n_user() < 0)
continue;
std::vector<cxxExchComp *> comps = exchange_ptr->Vectorize();
// check elements
for (size_t j = 0; j < comps.size(); j++)
{
if (comps[j]->Get_rate_name().size() == 0)
continue;
/* First find exchange master species */
cxxNameDouble nd = comps[j]->Get_totals();
cxxNameDouble::iterator kit = nd.begin();
bool found_exchange = false;
for (; kit != nd.end(); kit++)
{
/* Find master species */
struct element *elt_ptr = element_store(kit->first.c_str());
if (elt_ptr == NULL || elt_ptr->master == NULL)
{
input_error++;
error_string = sformatf( "Master species not in data "
"base for %s, skipping element.",
elt_ptr->name);
error_msg(error_string, CONTINUE);
continue;
}
if (elt_ptr->master->type == EX)
found_exchange = true;;
}
if (!found_exchange)
{
input_error++;
error_string = sformatf(
"Exchange formula does not contain an exchange master species, %s",
comps[j]->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
/* Now find associated kinetic reaction ... */
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, exchange_ptr->Get_n_user())) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use exchange related to kinetic reaction, %s",
exchange_ptr->Get_n_user(), comps[j]->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
size_t k;
for (k = 0; k < kinetics_ptr->Get_kinetics_comps().size(); k++)
{
if (strcmp_nocase
(comps[j]->Get_rate_name().c_str(),
kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to exchanger, %s, not found in KINETICS %d",
comps[j]->Get_rate_name().c_str(), comps[j]->Get_formula().c_str(), exchange_ptr->Get_n_user());
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
comps[j]->Set_rate_name(kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str());
/* make exchanger concentration proportional to mineral ... */
conc = kinetics_ptr->Get_kinetics_comps()[k].Get_m() * comps[j]->Get_phase_proportion();
count_elts = 0;
paren_count = 0;
{
char * temp_formula = string_duplicate(comps[j]->Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comps[j]->Set_totals(elt_list_NameDouble());
/*
* No check on availability of exchange elements
*/
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_kin_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion
@ -3014,6 +3120,7 @@ tidy_kin_exchange(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_min_exchange(void)
@ -3374,7 +3481,239 @@ tidy_min_surface(void)
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_kin_surface(void)
/* ---------------------------------------------------------------------- */
/*
* If surface is related to mineral, surface amount is
* set in proportion
*/
{
int i, j, n, jj, l;
cxxKinetics *kinetics_ptr;
struct surface_comp *comp_ptr;
struct master *master_ptr;
struct phase *phase_ptr;
//char token[MAX_LENGTH];
char *ptr;
LDBLE conc;
struct elt_list *elt_list_kinetics;
int count_elts_kinetics;
n = -999;
comp_ptr = NULL;
for (i = 0; i < count_surface; i++)
{
if (surface[i].new_def == FALSE)
continue;
if (surface[i].n_user < 0)
continue;
for (j = 0; j < surface[i].count_comps; j++)
{
if (surface[i].comps[j].rate_name == NULL)
continue;
comp_ptr = &surface[i].comps[j];
comp_ptr->master = NULL;
n = surface[i].n_user;
/* First find surface master species */
int k;
for (k = 0; comp_ptr->totals[k].elt != NULL; k++)
{
/* Find master species */
master_ptr = comp_ptr->totals[k].elt->master;
if (master_ptr == NULL)
{
input_error++;
error_string = sformatf( "Master species not in data "
"base for %s, skipping element.",
comp_ptr->totals[k].elt->name);
error_msg(error_string, CONTINUE);
continue;
}
if (master_ptr->type != SURF)
continue;
comp_ptr->master = master_ptr;
break;
}
if (comp_ptr->master == NULL)
{
input_error++;
error_string = sformatf(
"Surface formula does not contain a surface master species, %s",
comp_ptr->formula);
error_msg(error_string, CONTINUE);
continue;
}
/* Now find the kinetic reaction on which surface depends... */
//if ((kinetics_ptr = kinetics_bsearch(n, &k)) == NULL)
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n)) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use surface related to kinetic reaction, %s",
n, comp_ptr->formula);
error_msg(error_string, CONTINUE);
continue;
}
for (k = 0; k < (int) kinetics_ptr->Get_kinetics_comps().size(); k++)
{
cxxKineticsComp *kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
if (strcmp_nocase
(comp_ptr->rate_name,
kin_comp_ptr->Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == (int) kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to surface, %s, not found in Kinetics %d",
comp_ptr->rate_name, comp_ptr->formula, n);
error_msg(error_string, CONTINUE);
continue;
}
cxxKineticsComp *kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
/* use database name for phase */
comp_ptr->rate_name = string_hsave(kin_comp_ptr->Get_rate_name().c_str());
/* make surface concentration proportional to mineral ... */
conc = kin_comp_ptr->Get_m() * comp_ptr->phase_proportion;
/* if (conc < MIN_RELATED_SURFACE) conc = 0.0; */
{
char * temp_formula = string_duplicate(comp_ptr->formula);
ptr = temp_formula;
count_elts = 0;
paren_count = 0;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comp_ptr->totals =
(struct elt_list *) free_check_null(comp_ptr->totals);
comp_ptr->totals = elt_list_save();
/* area */
surface[i].charge[comp_ptr->charge].grams =
kin_comp_ptr->Get_m();
}
/*
* check on elements
*/
/* Go through each kinetic reaction, add all related surface compositions
* check for negative values
*/
if (surface[i].related_rate == FALSE)
continue;
//kinetics_ptr = kinetics_bsearch(n, &k);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n);
for (size_t k = 0; k < kinetics_ptr->Get_kinetics_comps().size(); k++)
{
cxxKineticsComp *kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
count_elts = 0;
paren_count = 0;
/* added in kinetics formula */
cxxNameDouble::iterator jit = kin_comp_ptr->Get_namecoef().begin();
//for (j = 0; j < kinetics_ptr->comps[k].count_list; j++)
for (; jit != kin_comp_ptr->Get_namecoef().end(); jit++)
{
std::string name = jit->first;
LDBLE coef = jit->second;
phase_ptr = NULL;
phase_ptr = phase_bsearch(name.c_str(), &jj, FALSE);
if (phase_ptr != NULL)
{
add_elt_list(phase_ptr->next_elt, 1.0);
}
else
{
char * temp_name = string_duplicate(name.c_str());
ptr = temp_name;
get_elts_in_species(&ptr, coef);
free_check_null(temp_name);
}
}
/* save kinetics formula */
if (count_elts > 0)
{
qsort(elt_list, (size_t) count_elts,
(size_t) sizeof(struct elt_list), elt_list_compare);
elt_list_combine();
}
elt_list_kinetics = elt_list_save();
count_elts_kinetics = count_elts;
/* get surface formulas */
count_elts = 0;
paren_count = 0;
for (j = 0; j < surface[i].count_comps; j++)
{
comp_ptr = &surface[i].comps[j];
if (comp_ptr->rate_name == NULL)
continue;
if (strcmp_nocase
(comp_ptr->rate_name,
kin_comp_ptr->Get_rate_name().c_str()) == 0)
{
char * temp_formula = string_duplicate( comp_ptr->formula);
ptr = temp_formula;
get_elts_in_species(&ptr,
-1 * comp_ptr->phase_proportion);
free_check_null(temp_formula);
}
}
if (count_elts > 0)
{
qsort(elt_list, (size_t) count_elts,
(size_t) sizeof(struct elt_list), elt_list_compare);
elt_list_combine();
}
for (j = 0; j < count_elts; j++)
{
if (elt_list[j].elt->primary->s->type <= H2O)
{
for (l = 0; l < count_elts_kinetics; l++)
{
if (elt_list[j].elt == elt_list_kinetics[l].elt)
{
break;
}
}
if (l == count_elts_kinetics)
{
input_error++;
error_string = sformatf(
"Stoichiometry of surface, %s * %g mol sites/mol reactant,\n\tmust be a subset of the formula defined for the related reactant %s.\n\tElement %s is not present in reactant formula.",
comp_ptr->formula,
(double) comp_ptr->phase_proportion,
comp_ptr->rate_name, elt_list[j].elt->name);
error_msg(error_string, CONTINUE);
}
else if (fabs(elt_list[j].coef) >
fabs(elt_list_kinetics[l].coef))
{
input_error++;
error_string = sformatf(
"Stoichiometry of surface, %s * %g mol sites/mol reactant,\n\tmust be a subset of the formula defined for the related reactant %s.\n\tCoefficient of element %s in surface exceeds amount present in reactant formula.",
comp_ptr->formula,
(double) comp_ptr->phase_proportion,
comp_ptr->rate_name, elt_list[j].elt->name);
error_msg(error_string, CONTINUE);
}
}
}
elt_list_kinetics =
(struct elt_list *) free_check_null(elt_list_kinetics);
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_kin_surface(void)
@ -3601,6 +3940,7 @@ tidy_kin_surface(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
ss_prep(LDBLE t, cxxSS *ss_ptr, int print)

26
phreeqc/transport.cpp
View File

@ -5,7 +5,7 @@
#include "GasPhase.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
transport(void)
@ -513,7 +513,8 @@ transport(void)
rate_sim_time = rate_sim_time_start + kin_time;
/* halftime kinetics for resident water in first cell ... */
if (kinetics_bsearch(first_c, &i) != NULL && count_cells > 1)
//if (kinetics_bsearch(first_c, &i) != NULL && count_cells > 1)
if (Utilities::Rxn_find(Rxn_kinetics_map, first_c) != NULL && count_cells > 1)
{
cell_no = first_c;
kin_time = kin_time_save / 2;
@ -1024,7 +1025,7 @@ int Phreeqc::
mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction)
/* ---------------------------------------------------------------------- */
{
int j, n, k, l;
int j, n, k;
LDBLE t_imm;
struct solution *ptr_imm, *ptr_m;
/*
@ -1078,7 +1079,8 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction)
set_and_run_wrapper(i, STAG, FALSE, -2, 0.0);
if (multi_Dflag == TRUE)
fill_spec(cell_no);
use.Set_kinetics_ptr(kinetics_bsearch(i, &l));
//use.Set_kinetics_ptr(kinetics_bsearch(i, &l));
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, i));
if (use.Get_kinetics_ptr() != NULL)
{
use.Set_n_kinetics_user(i);
@ -1325,9 +1327,9 @@ int Phreeqc::
set_initial_moles(int i)
/* ---------------------------------------------------------------------- */
{
struct kinetics *kinetics_ptr;
cxxKinetics *kinetics_ptr;
char token[MAX_LENGTH], token1[MAX_LENGTH], *ptr;
int j, k, l, n;
int j, k, l;
/*
* Pure phase assemblage
*/
@ -1363,11 +1365,17 @@ set_initial_moles(int i)
/*
* Kinetics
*/
kinetics_ptr = kinetics_bsearch(i, &n);
//kinetics_ptr = kinetics_bsearch(i, &n);
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, i);
if (kinetics_ptr != NULL)
{
for (j = 0; j < kinetics_ptr->count_comps; j++)
kinetics_ptr->comps[j].initial_moles = kinetics_ptr->comps[j].m;
//for (j = 0; j < kinetics_ptr->count_comps; j++)
// kinetics_ptr->comps[j].initial_moles = kinetics_ptr->comps[j].m;
for (j = 0; j < (int) kinetics_ptr->Get_kinetics_comps().size(); j++)
{
cxxKineticsComp *kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
kinetics_comp_ptr->Set_initial_moles(kinetics_comp_ptr->Get_m());
}
}
/*
* Solid solutions