From 694c07488e9acaa7cdb7de383d4aa5d0c56a38c6 Mon Sep 17 00:00:00 2001 From: David L Parkhurst Date: Thu, 14 Jan 2010 16:26:14 +0000 Subject: [PATCH] Bugs with kinetics in PHAST. The formula was not correct when converted to a phreeqcpp class, only the last instance of an element was kept and the rest were lost. Modified NameDouble to accumulate the sum of coefficients for each element when converting from a PHREEQC structure. The order of the kinetics components was not maintained when converted to a phreeqcpp class. Components were in alphabetical order. This posed a problem with get and put statements because the gets could end up before the puts. Changed from map of components to list of components, which should maintain the order correctly. Problems with modify solution and redox elements. Laurin pointed out that adding N may leave all the N(x) in place, which increases the total N. Revised read_solution to merge the valence states. If N defined, then all N(x) are removed. If N() defined, then N is removed. Still need to test. git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@3927 1feff8c3-07ed-0310-ac33-dd36852eb9cd --- KineticsComp.cxx | 28 +++++++-------- KineticsComp.h | 2 +- NameDouble.cxx | 88 ++++++++++++++++++++++++++++++++++++++++++++++++ NameDouble.h | 3 ++ Solution.cxx | 26 ++++++++++---- cxxKinetics.cxx | 66 +++++++++++++++++++++++++++++------- cxxKinetics.h | 3 +- 7 files changed, 182 insertions(+), 34 deletions(-) diff --git a/KineticsComp.cxx b/KineticsComp.cxx index f777a0c9..0d015b8b 100644 --- a/KineticsComp.cxx +++ b/KineticsComp.cxx @@ -59,7 +59,7 @@ cxxKineticsComp::~cxxKineticsComp() } struct kinetics_comp * - cxxKineticsComp::cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::map < std::string, cxxKineticsComp > + cxxKineticsComp::cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::list < cxxKineticsComp > &el) // // Builds kinetics_comp structure from of cxxKineticsComp @@ -72,20 +72,20 @@ struct kinetics_comp * P_INSTANCE_POINTER malloc_error(); int i = 0; - for (std::map < std::string, cxxKineticsComp >::iterator it = el.begin(); + for (std::list < cxxKineticsComp >::iterator it = el.begin(); it != el.end(); ++it) { - if ((*it).second.rate_name.size() == 0) + if ((*it).rate_name.size() == 0) kinetics_comp_ptr[i].rate_name = NULL; else - kinetics_comp_ptr[i].rate_name = P_INSTANCE_POINTER string_hsave((*it).second.rate_name.c_str()); - kinetics_comp_ptr[i].list = (*it).second.namecoef.name_coef(P_INSTANCE); - kinetics_comp_ptr[i].count_list = (int) (*it).second.namecoef.size(); - kinetics_comp_ptr[i].tol = (*it).second.tol; - kinetics_comp_ptr[i].m = (*it).second.m; + kinetics_comp_ptr[i].rate_name = P_INSTANCE_POINTER string_hsave((*it).rate_name.c_str()); + kinetics_comp_ptr[i].list = (*it).namecoef.name_coef(P_INSTANCE); + kinetics_comp_ptr[i].count_list = (int) (*it).namecoef.size(); + kinetics_comp_ptr[i].tol = (*it).tol; + kinetics_comp_ptr[i].m = (*it).m; kinetics_comp_ptr[i].initial_moles = 0.; - kinetics_comp_ptr[i].m0 = (*it).second.m0; - kinetics_comp_ptr[i].moles = (*it).second.moles; + kinetics_comp_ptr[i].m0 = (*it).m0; + kinetics_comp_ptr[i].moles = (*it).moles; kinetics_comp_ptr[i].count_c_params = 0; kinetics_comp_ptr[i].c_params = NULL; /* @@ -93,15 +93,15 @@ struct kinetics_comp * kinetics_comp_ptr[i].d_params = NULL; */ - kinetics_comp_ptr[i].count_d_params = (int) (*it).second.d_params.size(); + kinetics_comp_ptr[i].count_d_params = (int) (*it).d_params.size(); kinetics_comp_ptr[i].d_params = NULL; - if ((*it).second.d_params.size() > 0) + if ((*it).d_params.size() > 0) { kinetics_comp_ptr[i].d_params = (double *) - P_INSTANCE_POINTER PHRQ_malloc((size_t) ((*it).second.d_params.size() * sizeof(double))); + P_INSTANCE_POINTER PHRQ_malloc((size_t) ((*it).d_params.size() * sizeof(double))); if (kinetics_comp_ptr[i].d_params == NULL) P_INSTANCE_POINTER malloc_error(); - std::copy((*it).second.d_params.begin(), (*it).second.d_params.end(), + std::copy((*it).d_params.begin(), (*it).d_params.end(), kinetics_comp_ptr[i].d_params); } i++; diff --git a/KineticsComp.h b/KineticsComp.h index 754034d7..acecc8c7 100644 --- a/KineticsComp.h +++ b/KineticsComp.h @@ -17,7 +17,7 @@ public: cxxKineticsComp(struct kinetics_comp *); ~cxxKineticsComp(); - static struct kinetics_comp *cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::map < std::string, cxxKineticsComp > &el); + static struct kinetics_comp *cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::list < cxxKineticsComp > &el); void dump_xml(std::ostream & os, unsigned int indent = 0) const; diff --git a/NameDouble.cxx b/NameDouble.cxx index c36fead0..2c16de98 100644 --- a/NameDouble.cxx +++ b/NameDouble.cxx @@ -114,6 +114,7 @@ cxxNameDouble::cxxNameDouble(struct master_activity *ma, int count, } this->type = type; } +/* cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count) // // constructor for cxxNameDouble from list of elt_list @@ -128,6 +129,30 @@ cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count) } this->type = ND_NAME_COEF; } +*/ +cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count) + // + // constructor for cxxNameDouble from list of elt_list + // +{ + int i; + for (i = 0; i < count; i++) + { + if (nc[i].name == NULL) + continue; + + if ((*this).find(nc[i].name) == (*this).end()) + { + (*this)[nc[i].name] = nc[i].coef; + } + else + { + (*this)[nc[i].name] = (*this).find(nc[i].name)->second + nc[i].coef; + } + + } + this->type = ND_NAME_COEF; +} cxxNameDouble::~cxxNameDouble() { @@ -462,7 +487,70 @@ cxxNameDouble::multiply(double extensive) it->second *= extensive; } } +void +cxxNameDouble::merge_redox(const cxxNameDouble & source) +// +// Merges source into this +// Accounts for possible conflicts between redox state and +// totals +// +{ + for (cxxNameDouble::const_iterator sit = source.begin(); sit != source.end(); sit++) + { + + std::string redox_name = sit->first; + std::string elt_name; + size_t pos = redox_name.find("("); + bool redox; + if (pos != std::string::npos) + { + redox = true; + elt_name = redox_name.substr(0, pos - 1); + } + else + { + redox = false; + elt_name = redox_name; + } + if (redox) + { + // Remove elt_name, if present + if ((*this).find(elt_name) != (*this).end()) + { + (*this).erase((*this).find(elt_name)); + } + // Put in redox name + (*this)[redox_name] = sit->second; + } + else + { + std::string substring; + substring.append(elt_name); + substring.append("("); + + // Remove all redox + bool deleted = true; + while (deleted) + { + deleted = false; + cxxNameDouble::iterator current = (*this).begin(); + for ( ; current != (*this).end(); current++) + { + if (current->first.find(substring) == 0) + { + (*this).erase(current); + deleted = true; + break; + } + } + } + // Put in elt name + (*this)[elt_name] = sit->second; + } + + } +} #ifdef USE_MPI void cxxNameDouble::mpi_pack(std::vector < int >&ints, diff --git a/NameDouble.h b/NameDouble.h index 09e0f541..e94c533e 100644 --- a/NameDouble.h +++ b/NameDouble.h @@ -66,6 +66,9 @@ class cxxNameDouble:public void multiply(double factor); + void + merge_redox(const cxxNameDouble & source); + void insert(char *str, double d) { diff --git a/Solution.cxx b/Solution.cxx index 0eecc1c3..f5a15ba5 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -630,13 +630,27 @@ cxxSolution::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check) continue; case 0: // totals - if (this->totals.read_raw(P_INSTANCE_COMMA parser, next_char) != - CParser::PARSER_OK) + //if (this->totals.read_raw(P_INSTANCE_COMMA parser, next_char) != + // CParser::PARSER_OK) + //{ + // parser.incr_input_error(); + // parser. + // error_msg("Expected element name and moles for totals.", + // CParser::OT_CONTINUE); + //} { - parser.incr_input_error(); - parser. - error_msg("Expected element name and moles for totals.", - CParser::OT_CONTINUE); + cxxNameDouble temp_totals; + if (temp_totals.read_raw(P_INSTANCE_COMMA parser, next_char) != CParser::PARSER_OK) + { + parser.incr_input_error(); + parser. + error_msg("Expected element name and moles for totals.", + CParser::OT_CONTINUE); + } + else + { + this->totals.merge_redox(temp_totals); + } } opt_save = 0; break; diff --git a/cxxKinetics.cxx b/cxxKinetics.cxx index 02d50055..9856cd8b 100644 --- a/cxxKinetics.cxx +++ b/cxxKinetics.cxx @@ -66,7 +66,8 @@ totals(kinetics_ptr->totals) { cxxKineticsComp ec(&(kinetics_ptr->comps[i])); std::string str(ec.get_rate_name()); - this->kineticsComps[str] = ec; + //this->kineticsComps[str] = ec; + this->kineticsComps.push_back(ec); } // steps @@ -231,12 +232,12 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent) const s_oss << "-cvode_order " << this->cvode_order << std::endl; // kineticsComps structures - for (std::map < std::string, cxxKineticsComp >::const_iterator it = + for (std::list < cxxKineticsComp >::const_iterator it = kineticsComps.begin(); it != kineticsComps.end(); ++it) { s_oss << indent1; s_oss << "-component" << std::endl; - (*it).second.dump_raw(s_oss, indent + 2); + (*it).dump_raw(s_oss, indent + 2); } // totals @@ -421,6 +422,28 @@ cxxKinetics::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check) reread.set_echo_file(CParser::EO_NONE); reread.set_echo_stream(CParser::EO_NONE); + std::list < cxxKineticsComp >::iterator kit; + bool found = false; + for (kit = this->kineticsComps.begin(); kit != this->kineticsComps.end(); kit++) + { + if (kit->get_rate_name() == ec.get_rate_name()) + { + found = true; + break; + } + } + if (found) + { + kit->read_raw(P_INSTANCE_COMMA reread, false); + } + else + { + cxxKineticsComp ec1; + ec1.read_raw(P_INSTANCE_COMMA reread, false); + std::string str(ec1.get_rate_name()); + this->kineticsComps.push_back(ec1); + } + /* if (this->kineticsComps.find(ec.get_rate_name()) != this->kineticsComps.end()) { cxxKineticsComp & comp = this->kineticsComps.find(ec.get_rate_name())->second; @@ -433,6 +456,7 @@ cxxKinetics::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check) std::string str(ec1.get_rate_name()); this->kineticsComps[str] = ec1; } + */ } useLastLine = true; break; @@ -554,11 +578,18 @@ cxxKinetics::mpi_pack(std::vector < int >&ints, { ints.push_back(this->n_user); ints.push_back((int) this->kineticsComps.size()); + std::list < cxxKineticsComp >::iterator it; + for (it = this->kineticsComps.begin(); it != this->kineticsComps.end(); it++) + { + (*it).mpi_pack(ints, doubles); + } + /* for (std::map < std::string, cxxKineticsComp >::iterator it = this->kineticsComps.begin(); it != this->kineticsComps.end(); it++) { (*it).second.mpi_pack(ints, doubles); } + */ ints.push_back((int) this->steps.size()); for (std::vector < double >::iterator it = this->steps.begin(); it != this->steps.end(); it++) @@ -588,8 +619,9 @@ cxxKinetics::mpi_unpack(int *ints, int *ii, double *doubles, int *dd) { cxxKineticsComp kc; kc.mpi_unpack(ints, &i, doubles, &d); - std::string str(kc.get_rate_name()); - this->kineticsComps[str] = kc; + //std::string str(kc.get_rate_name()); + //this->kineticsComps[str] = kc; + this->kineticsComps.push_back(kc); } n = ints[i++]; this->steps.clear(); @@ -664,21 +696,31 @@ cxxKinetics::add(const cxxKinetics & addee, double extensive) if (extensive == 0.0) return; //std::map < std::string, cxxKineticsComp> kineticsComps; - for (std::map < std::string, cxxKineticsComp >::const_iterator itadd = + for (std::list < cxxKineticsComp >::const_iterator itadd = addee.kineticsComps.begin(); itadd != addee.kineticsComps.end(); ++itadd) { - std::map < std::string, cxxKineticsComp >::iterator it = this->kineticsComps.find((*itadd).first); - if (it != this->kineticsComps.end()) + bool found(false); + std::list < cxxKineticsComp >::iterator it; + for (it = this->kineticsComps.begin(); it != this->kineticsComps.end(); it++) { - (*it).second.add((*itadd).second, extensive); + if ((*it).get_rate_name() == (*itadd).get_rate_name()) + { + found = true; + break; + } + } + + //std::map < std::string, cxxKineticsComp >::iterator it = this->kineticsComps.find((*itadd).first); + if (found) + { + (*it).add((*itadd), extensive); } else { - cxxKineticsComp entity = (*itadd).second; + cxxKineticsComp entity = (*itadd); entity.multiply(extensive); - std::string str(entity.get_rate_name()); - this->kineticsComps[str] = entity; + this->kineticsComps.push_back(entity); } } //std::vector steps; diff --git a/cxxKinetics.h b/cxxKinetics.h index 750b4d03..014b911d 100644 --- a/cxxKinetics.h +++ b/cxxKinetics.h @@ -43,7 +43,8 @@ class cxxKinetics:public cxxNumKeyword void add(const cxxKinetics & addee, double extensive); protected: - std::map < std::string, cxxKineticsComp > kineticsComps; + //std::map < std::string, cxxKineticsComp > kineticsComps; + std::list < cxxKineticsComp > kineticsComps; std::vector < double >steps; cxxNameDouble totals; double step_divide;