From 2b442eef3006739d8bcee7d0e1d4992403c1b41e Mon Sep 17 00:00:00 2001 From: David L Parkhurst Date: Tue, 30 Jan 2007 01:34:18 +0000 Subject: [PATCH] git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@1723 1feff8c3-07ed-0310-ac33-dd36852eb9cd --- ExchComp.cxx | 2 +- Exchange.h | 14 ++-- NameDouble.cxx | 68 +++------------- NameDouble.h | 4 +- Solution.cxx | 172 ++++++---------------------------------- Solution.h | 22 +++-- SolutionIsotopeList.cxx | 93 +++++++++++----------- StorageBin.cxx | 4 +- StorageBin.h | 14 +--- 9 files changed, 107 insertions(+), 286 deletions(-) diff --git a/ExchComp.cxx b/ExchComp.cxx index f6373344..d7848d0f 100644 --- a/ExchComp.cxx +++ b/ExchComp.cxx @@ -102,7 +102,7 @@ cxxExchComp::cxxExchComp(std::vector &ec_vector, std::vectorla += it_ec->la*intensive; this->charge_balance += it_ec->charge_balance*extensive; this->phase_proportion += it_ec->phase_proportion*intensive; - this->totals.add_extensive(it_ec->totals, extensive); + this->totals.add(it_ec->totals, extensive); it_ec++; it_f++; } diff --git a/Exchange.h b/Exchange.h index 22d2f0fa..ac7c1a8c 100644 --- a/Exchange.h +++ b/Exchange.h @@ -1,15 +1,17 @@ #if !defined(EXCHANGE_H_INCLUDED) #define EXCHANGE_H_INCLUDED +#include "NumKeyword.h" #define EXTERNAL extern #include "global.h" #include // assert #include // std::map #include // std::string -#include // std::list +#include // std::list +#include // std::vector +#include "char_star.h" #include "ExchComp.h" -#include "cxxMix.h" class cxxExchange : public cxxNumKeyword { @@ -17,13 +19,14 @@ class cxxExchange : public cxxNumKeyword public: cxxExchange(); cxxExchange(struct exchange *); - cxxExchange(const std::map &exchange_map, cxxMix &mx); ~cxxExchange(); struct exchange *cxxExchange2exchange(); struct exch_comp *cxxExchComp2exch_comp(); + void dump_xml(std::ostream& os, unsigned int indent = 0)const; + void dump_raw(std::ostream& s_oss, unsigned int indent)const; void read_raw(CParser& parser); @@ -47,11 +50,6 @@ public: void mpi_pack(std::vector& ints, std::vector& doubles); void mpi_unpack(int *ints, int *ii, double *doubles, int *dd); #endif -private: - void zero(); - void add(const cxxExchange &ex, double intensive, double extensive); - // not written - void dump_xml(std::ostream& os, unsigned int indent = 0)const; protected: std::list exchComps; diff --git a/NameDouble.cxx b/NameDouble.cxx index e63b47a9..5561713d 100644 --- a/NameDouble.cxx +++ b/NameDouble.cxx @@ -275,63 +275,21 @@ CParser::STATUS_TYPE cxxNameDouble::read_raw(CParser& parser, std::istream::pos_ return CParser::PARSER_OK; } -void cxxNameDouble::add_extensive(const cxxNameDouble &addee, double factor) -// -// Sums two name doubles, this + factor*nd2 -// +void cxxNameDouble::add(const cxxNameDouble &old, double factor) + // + // constructor for cxxNameDouble from list of elt_list + // { - if (factor == 0) return; - assert (factor > 0.0); - for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++) - { - cxxNameDouble::iterator current = (*this).find(it->first); - if (current != (*this).end()) - { - (*this)[it->first] = current->second + it->second * factor; - } else { - (*this)[it->first] = it->second * factor; - } - } -} -void cxxNameDouble::add_intensive(const cxxNameDouble &addee, double f1, double f2) -// -// Sums two name doubles, this*f1 + f2*nd2 -// -{ - assert(f1 > 0 && f2 > 0); - for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++) - { - cxxNameDouble::iterator current = (*this).find(it->first); - if (current != (*this).end()) - { - (*this)[it->first] = f1*current->second + f2*it->second; - } else { - (*this)[it->first] = f2 * it->second; - } - } -} -void cxxNameDouble::add_log_activities(const cxxNameDouble &addee, double f1, double f2) -// -// Sums two name doubles, this*f1 + f2*nd2, assuming log values -// -{ - assert (f1 > 0 && f2 > 0); - for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++) - { - cxxNameDouble::iterator current = (*this).find(it->first); - if (current != (*this).end()) - { - double a1 = pow(10., current->second); - double a2 = pow(10., it->second); - (*this)[it->first] = log10(f1*a1 + f2*a2); - } else { - //double a2 = pow(10. it->second); - //(*this)[it->first] = log10(f2 * a2); - (*this)[it->first] = it->second + log10(f2); - } - } -} + for (cxxNameDouble::const_iterator it = old.begin(); it != old.end(); it++) { + cxxNameDouble::iterator current = (*this).find(it->first); + if (current != (*this).end()) { + (*this)[it->first] = current->second + it->second * factor; + } else { + (*this)[it->first] = it->second * factor; + } + } +} void cxxNameDouble::add(char * key, double total) // // add to total for a specified element diff --git a/NameDouble.h b/NameDouble.h index 73e74584..7623c2f5 100644 --- a/NameDouble.h +++ b/NameDouble.h @@ -44,9 +44,7 @@ public: CParser::STATUS_TYPE read_raw(CParser& parser, std::istream::pos_type& pos); - void add_extensive(const cxxNameDouble &old, double factor); - void add_intensive(const cxxNameDouble &addee, double fthis, double f2); - void add_log_activities(const cxxNameDouble &addee, double fthis, double f2); + void add(const cxxNameDouble &old, double factor); void add(char * key, double total); void insert(char *str, double d) { diff --git a/Solution.cxx b/Solution.cxx index b697ae36..cb1e65c7 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -27,21 +27,20 @@ cxxSolution::cxxSolution() // : cxxNumKeyword() { - this->tc = 25.0; - this->ph = 7.0; - this->pe = 4.0; - this->mu = 1e-7; - this->ah2o = 1.0; - this->total_h = 111.1; - this->total_o = 55.55; - this->cb = 0.0; - this->mass_water = 1.0; - this->total_alkalinity = 0.0; - this->totals.type = cxxNameDouble::ND_ELT_MOLES; - this->master_activity.type = cxxNameDouble::ND_SPECIES_LA; - this->species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; + tc = 25.0; + ph = 7.0; + pe = 4.0; + mu = 1e-7; + ah2o = 1.0; + total_h = 111.1; + total_o = 55.55; + cb = 0.0; + mass_water = 1.0; + total_alkalinity = 0.0; + totals.type = cxxNameDouble::ND_ELT_MOLES; + master_activity.type = cxxNameDouble::ND_SPECIES_LA; + species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; } -/* cxxSolution::cxxSolution(double zero) // // empty cxxSolution constructor @@ -63,7 +62,6 @@ cxxSolution::cxxSolution(double zero) master_activity.type = cxxNameDouble::ND_SPECIES_LA; species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; } -*/ cxxSolution::cxxSolution(struct solution *solution_ptr) // // constructor for cxxSolution from struct solution @@ -77,18 +75,18 @@ isotopes(solution_ptr) { this->set_description(solution_ptr->description); - this->n_user = solution_ptr->n_user; - this->n_user_end = solution_ptr->n_user_end; - this->tc = solution_ptr->tc; - this->ph = solution_ptr->ph; - this->pe = solution_ptr->solution_pe; - this->mu = solution_ptr->mu; - this->ah2o = solution_ptr->ah2o; - this->total_h = solution_ptr->total_h; - this->total_o = solution_ptr->total_o; - this->cb = solution_ptr->cb; - this->mass_water = solution_ptr->mass_water; - this->total_alkalinity = solution_ptr->total_alkalinity; + n_user = solution_ptr->n_user; + n_user_end = solution_ptr->n_user_end; + tc = solution_ptr->tc; + ph = solution_ptr->ph; + pe = solution_ptr->solution_pe; + mu = solution_ptr->mu; + ah2o = solution_ptr->ah2o; + total_h = solution_ptr->total_h; + total_o = solution_ptr->total_o; + cb = solution_ptr->cb; + mass_water = solution_ptr->mass_water; + total_alkalinity = solution_ptr->total_alkalinity; // Totals filled in constructor, just save description and moles @@ -103,83 +101,6 @@ isotopes(solution_ptr) // Species_gamma in constructor } -#ifdef SKIP -cxxSolution::cxxSolution( const std::map& solutions, cxxMix &mix) - // - // constructor for cxxSolution from struct solution - //cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap) - -{ -/* -* mixes solutions based on cxxMix structure, returns new solution -* return solution must be freed by calling method -*/ - double intensive, extensive; -/* -* Zero out global solution data -*/ - this->zero(); -/* -* Determine sum of mixing fractions -*/ - extensive = 0.0; - - std::map *mixcomps = mix.comps(); - - std::map::const_iterator it; - for (it = mixcomps->begin(); it != mixcomps->end(); it++) { - extensive += it->second; - } -/* -* Add solutions -*/ - for (it = mixcomps->begin(); it != mixcomps->end(); it++) - { - const cxxSolution *cxxsoln_ptr1 = &(solutions.find(it->first)->second); - if (cxxsoln_ptr1 == NULL) - { - sprintf(error_string, "Solution %d not found in mix_cxxSolutions.", it->first); - error_msg(error_string, CONTINUE); - input_error++; - return; - } - intensive = it->second/extensive; - double ext = it->second; - this->add(*cxxsoln_ptr1, intensive, it->second); - } -} -#endif - -cxxSolution::cxxSolution( const std::map& solutions, cxxMix &mix) -// -// constructor for cxxSolution from mixture of solutions -// - -{ - -// -// Zero out solution data -// - this->zero(); -// -// Mix solutions -// - std::map *mixcomps = mix.comps(); - std::map::const_iterator it; - for (it = mixcomps->begin(); it != mixcomps->end(); it++) - { - const cxxSolution *cxxsoln_ptr1 = &(solutions.find(it->first)->second); - if (cxxsoln_ptr1 == NULL) - { - sprintf(error_string, "Solution %d not found in mix_cxxSolutions.", it->first); - error_msg(error_string, CONTINUE); - input_error++; - return; - } - this->add(*cxxsoln_ptr1, it->second); - } -} - #ifdef SKIP cxxSolution::cxxSolution(const cxxSolution &old, double intensive, double extensive) // @@ -748,23 +669,6 @@ void cxxSolution::read_raw(CParser& parser) return; } -void cxxSolution::zero() -{ - this->tc = 0.0; - this->ph = 0.0; - this->pe = 0.0; - this->mu = 0.0; - this->ah2o = 0.0; - this->total_h = 0.0; - this->total_o = 0.0; - this->cb = 0.0; - this->mass_water = 0.0; - this->total_alkalinity = 0.0; - this->totals.type = cxxNameDouble::ND_ELT_MOLES; - this->master_activity.type = cxxNameDouble::ND_SPECIES_LA; - this->species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; -} -#ifdef SKIP void cxxSolution::add(const cxxSolution &old, double intensive, double extensive) // // Add existing solution to "this" solution @@ -818,31 +722,7 @@ void cxxSolution::add(const cxxSolution &old, double intensive, double extensive cxxNameDouble species_gamma; */ } -#endif -void cxxSolution::add(const cxxSolution &addee, double extensive) - // - // Add existing solution to "this" solution - // -{ - double ext1 = this->mass_water; - double ext2 = addee.mass_water * extensive; - double f1 = ext1/(ext1 + ext2); - double f2 = ext2/(ext1 + ext2); - this->tc = f1*this->tc + f2*addee.tc; - this->ph = f1*this->ph + f2*addee.ph; - this->pe = f1*this->pe + f2*addee.pe; - this->mu = f1*this->mu + f2*addee.mu; - this->ah2o = f1*this->mu + f2*addee.ah2o; - this->total_h += addee.total_h * extensive; - this->total_o += addee.total_o * extensive; - this->cb += addee.cb * extensive; - this->mass_water += addee.mass_water * extensive; - this->total_alkalinity += addee.total_alkalinity * extensive; - this->totals.add_extensive(addee.totals, extensive); - this->master_activity.add_log_activities(addee.master_activity, f1, f2); - this->species_gamma.add_intensive(addee.species_gamma, f1, f2); - this->isotopes.add(addee.isotopes, f2, extensive); -} + double cxxSolution::get_total(char *string)const { cxxNameDouble::const_iterator it = this->totals.find(string); diff --git a/Solution.h b/Solution.h index 23f01edd..52c3983c 100644 --- a/Solution.h +++ b/Solution.h @@ -11,7 +11,7 @@ #include // assert #include // std::map #include // std::string - +#include // std::list #include // std::vector #include "char_star.h" @@ -21,11 +21,12 @@ class cxxSolution : public cxxNumKeyword public: cxxSolution(); - //cxxSolution(double zero); + cxxSolution(double zero); cxxSolution(struct solution *); cxxSolution(const std::map &solution_map, cxxMix &mx); - //cxxSolution(const cxxSolution &old, double intensive, double extensive); - //cxxSolution(const cxxSolution&); + cxxSolution(const cxxSolution &old, double intensive, double extensive); + + //cxxSolution(const cxxSolution&); ~cxxSolution(); //static cxxSolution& read(CParser& parser); @@ -76,10 +77,13 @@ public: struct solution *cxxSolution2solution(); - void dump_raw(std::ostream& s_oss, unsigned int indent)const; + void dump_xml(std::ostream& os, unsigned int indent = 0)const; + + void dump_raw(std::ostream& s_oss, unsigned int indent)const; void read_raw(CParser& parser); + void add(const cxxSolution &sol, double intensive, double extensive); #ifdef USE_MPI void mpi_pack(std::vector& ints, std::vector& doubles); @@ -87,14 +91,6 @@ public: void mpi_send(int task_number); void mpi_recv(int task_number); #endif - -private: - void zero(); - //void add(const cxxSolution &sol, double intensive, double extensive); - void add(const cxxSolution &addee, double extensive); - // not checked - void dump_xml(std::ostream& os, unsigned int indent = 0)const; - protected: double tc; double ph; diff --git a/SolutionIsotopeList.cxx b/SolutionIsotopeList.cxx index 0bc81085..88a0d135 100644 --- a/SolutionIsotopeList.cxx +++ b/SolutionIsotopeList.cxx @@ -32,55 +32,54 @@ cxxSolutionIsotopeList::cxxSolutionIsotopeList(struct solution *solution_ptr) } void cxxSolutionIsotopeList::add(cxxSolutionIsotopeList old, double intensive, double extensive) { - for (cxxSolutionIsotopeList::const_iterator itold = old.begin(); itold != old.end(); ++itold) { - //for (std::list ::const_iterator itold = old.isotopes.begin(); itold != old.isotopes.end(); ++itold) { - bool found = false; - for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) - { - //for (std::list ::iterator it = this->isotopes.begin(); it != this->isotopes.end(); ++it) { - if ((it->isotope_number == itold->isotope_number) && - (it->elt_name == itold->elt_name) && - (it->isotope_name == itold->isotope_name)) - { - it->total += itold->total * extensive; - it->ratio += itold->ratio * intensive; - it->ratio_uncertainty += itold->ratio_uncertainty * intensive; - it->ratio_uncertainty_defined = (it->ratio_uncertainty_defined || itold->ratio_uncertainty_defined); - found = true; - break; - } - } - if (!found) - { - cxxSolutionIsotope iso; - iso.total = itold->total * extensive; - iso.ratio = itold->ratio * intensive; - iso.ratio_uncertainty = itold->ratio_uncertainty * intensive; - iso.ratio_uncertainty_defined = itold->ratio_uncertainty_defined; - this->push_back(iso); - } - } + + for (cxxSolutionIsotopeList::const_iterator itold = old.begin(); itold != old.end(); ++itold) { + + //for (std::list ::const_iterator itold = old.isotopes.begin(); itold != old.isotopes.end(); ++itold) { + bool found = false; + for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) { + //for (std::list ::iterator it = this->isotopes.begin(); it != this->isotopes.end(); ++it) { + if ((it->isotope_number == itold->isotope_number) && + (it->elt_name == itold->elt_name) && + (it->isotope_name == itold->isotope_name)) { + it->total += itold->total * extensive; + it->ratio += itold->ratio * intensive; + it->ratio_uncertainty += itold->ratio_uncertainty * intensive; + it->ratio_uncertainty_defined = (it->ratio_uncertainty_defined || itold->ratio_uncertainty_defined); + found = true; + break; + } + } + if (!found) { + cxxSolutionIsotope iso; + iso.total = itold->total * extensive; + iso.ratio = itold->ratio * intensive; + iso.ratio_uncertainty = itold->ratio_uncertainty * intensive; + iso.ratio_uncertainty_defined = itold->ratio_uncertainty_defined; + this->push_back(iso); + } + } } struct isotope * cxxSolutionIsotopeList::cxxSolutionIsotopeList2isotope() { - struct isotope *iso; - if (this->size() <= 0) { - return NULL; - } else { - iso = (struct isotope *) PHRQ_malloc((size_t) ((this->size()) * sizeof(struct isotope))); - if (iso == NULL) malloc_error(); - int i = 0; - for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) { - iso[i].isotope_number = it->isotope_number; - iso[i].elt_name = it->elt_name; - iso[i].total = it->total; - iso[i].ratio = it->ratio; - iso[i].ratio_uncertainty = it->ratio_uncertainty; - iso[i].master = it->master(); - iso[i].primary = it->primary(); - i++; - } - } - return(iso); + struct isotope *iso; + if (this->size() <= 0) { + return NULL; + } else { + iso = (struct isotope *) PHRQ_malloc((size_t) ((this->size()) * sizeof(struct isotope))); + if (iso == NULL) malloc_error(); + int i = 0; + for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) { + iso[i].isotope_number = it->isotope_number; + iso[i].elt_name = it->elt_name; + iso[i].total = it->total; + iso[i].ratio = it->ratio; + iso[i].ratio_uncertainty = it->ratio_uncertainty; + iso[i].master = it->master(); + iso[i].primary = it->primary(); + i++; + } + } + return(iso); } diff --git a/StorageBin.cxx b/StorageBin.cxx index 6eeb6da7..d8d9f3b0 100644 --- a/StorageBin.cxx +++ b/StorageBin.cxx @@ -495,7 +495,7 @@ void cxxStorageBin::remove(int n) // Surface this->Surfaces.erase(n); } -#ifdef SKIP + cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap) { @@ -536,7 +536,7 @@ cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap) } return(cxxsoln_ptr); } -#endif + struct system *cxxStorageBin::cxxStorageBin2system(int n) // // make a system from storagebin diff --git a/StorageBin.h b/StorageBin.h index aa8806e0..3b19c292 100644 --- a/StorageBin.h +++ b/StorageBin.h @@ -42,8 +42,8 @@ public: } return (NULL); } - void setSolution(int n_user, cxxSolution &entity) { - Solutions[n_user] = entity; + void setSolution(int n_user, cxxSolution *entity) { + Solutions[n_user] = *entity; } void removeSolution(int n_user) { Solutions.erase(n_user); @@ -137,17 +137,9 @@ public: struct system *cxxStorageBin2system(int i); - //cxxSolution *mix_cxxSolutions(cxxMix &mixmap); + cxxSolution *mix_cxxSolutions(cxxMix &mixmap); cxxExchange *mix_cxxExchange(cxxMix &mixmap); - const std::map& getSolutions()const {return this->Solutions;}; - const std::map& getExchangers()const {return this->Exchangers;}; - const std::map& getGasPhases()const {return this->GasPhases;}; - const std::map& getKinetics()const {return this->Kinetics;}; - const std::map& getPPassemblages()const {return this->PPassemblages;}; - const std::map& getSSassemblages()const {return this->SSassemblages;}; - const std::map& getSurfaces()const {return this->Surfaces;}; - #ifdef USE_MPI void mpi_send(int n, int task_number); void mpi_recv(int task_number);