diff --git a/ExchComp.cxx b/ExchComp.cxx index d7848d0f..f6373344 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(it_ec->totals, extensive); + this->totals.add_extensive(it_ec->totals, extensive); it_ec++; it_f++; } diff --git a/Exchange.h b/Exchange.h index ac7c1a8c..22d2f0fa 100644 --- a/Exchange.h +++ b/Exchange.h @@ -1,17 +1,15 @@ #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::vector +#include // std::list -#include "char_star.h" #include "ExchComp.h" +#include "cxxMix.h" class cxxExchange : public cxxNumKeyword { @@ -19,14 +17,13 @@ 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); @@ -50,6 +47,11 @@ 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 5561713d..e63b47a9 100644 --- a/NameDouble.cxx +++ b/NameDouble.cxx @@ -275,21 +275,63 @@ CParser::STATUS_TYPE cxxNameDouble::read_raw(CParser& parser, std::istream::pos_ return CParser::PARSER_OK; } -void cxxNameDouble::add(const cxxNameDouble &old, double factor) - // - // constructor for cxxNameDouble from list of elt_list - // +void cxxNameDouble::add_extensive(const cxxNameDouble &addee, double factor) +// +// Sums two name doubles, this + factor*nd2 +// { - 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; - } - } - + 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); + } + } +} + void cxxNameDouble::add(char * key, double total) // // add to total for a specified element diff --git a/NameDouble.h b/NameDouble.h index 7623c2f5..73e74584 100644 --- a/NameDouble.h +++ b/NameDouble.h @@ -44,7 +44,9 @@ public: CParser::STATUS_TYPE read_raw(CParser& parser, std::istream::pos_type& pos); - void add(const cxxNameDouble &old, double factor); + 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(char * key, double total); void insert(char *str, double d) { diff --git a/Solution.cxx b/Solution.cxx index cb1e65c7..b697ae36 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -27,20 +27,21 @@ cxxSolution::cxxSolution() // : cxxNumKeyword() { - 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; + 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; } +/* cxxSolution::cxxSolution(double zero) // // empty cxxSolution constructor @@ -62,6 +63,7 @@ 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 @@ -75,18 +77,18 @@ isotopes(solution_ptr) { this->set_description(solution_ptr->description); - 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; + 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; // Totals filled in constructor, just save description and moles @@ -101,6 +103,83 @@ 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) // @@ -669,6 +748,23 @@ 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 @@ -722,7 +818,31 @@ 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 52c3983c..23f01edd 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,12 +21,11 @@ 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); @@ -77,13 +76,10 @@ public: struct solution *cxxSolution2solution(); - void dump_xml(std::ostream& os, unsigned int indent = 0)const; - - void dump_raw(std::ostream& s_oss, unsigned int indent)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); @@ -91,6 +87,14 @@ 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 88a0d135..0bc81085 100644 --- a/SolutionIsotopeList.cxx +++ b/SolutionIsotopeList.cxx @@ -32,54 +32,55 @@ 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 d8d9f3b0..6eeb6da7 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 3b19c292..aa8806e0 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,9 +137,17 @@ 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);