diff --git a/ExchComp.cxx b/ExchComp.cxx index 2b92367d..d51f3a27 100644 --- a/ExchComp.cxx +++ b/ExchComp.cxx @@ -25,13 +25,14 @@ cxxExchComp::cxxExchComp() // { moles = 0.0; + formula_totals.type = cxxNameDouble::ND_ELT_MOLES; + totals.type = cxxNameDouble::ND_ELT_MOLES; la = 0.0; charge_balance = 0.0; phase_name = NULL; phase_proportion = 0.0; rate_name = NULL; - totals.type = cxxNameDouble::ND_ELT_MOLES; - formula_totals.type = cxxNameDouble::ND_ELT_MOLES; + formula_z = 0.0; } cxxExchComp::cxxExchComp(struct exch_comp *exch_comp_ptr) @@ -53,6 +54,58 @@ totals(exch_comp_ptr->totals) rate_name = exch_comp_ptr->rate_name; formula_z = exch_comp_ptr->formula_z; } +cxxExchComp::cxxExchComp(std::vector &ec_vector, std::vector &f_vector) + // + // constructor for cxxExchComp from struct exch_comp + // +{ + if (ec_vector.size() <= 0) return; + // + // check consistency + // + std::vector::iterator it_f; + std::vector::iterator it_ec; + // set fixed variables + it_ec = ec_vector.begin(); + this->formula = it_ec->formula; + this->formula_totals = it_ec->formula_totals; + this->formula_z = it_ec->formula_z; + this->phase_name = it_ec->phase_name; + this->rate_name = it_ec->rate_name; + it_ec++; + for ( ; it_ec != ec_vector.end(); it_ec++) { + if (it_ec->formula != this->formula || + it_ec->formula_z != this->formula_z || + it_ec->phase_name != this->phase_name || + this->rate_name != this->rate_name) { + error_msg("Mixing exchange components. Formula, z, phase_name, or rate_name did not match", STOP); + } + } + // calculate sum of extensive factors + double sum_extensive = 0; + for (it_f = f_vector.begin(); it_f != f_vector.end(); it_f++) { + sum_extensive += *it_f; + } + this->moles = 0; + this->la = 0; + this->charge_balance = 0; + this->phase_proportion = 0; + this->totals.clear(); + this->totals.type = cxxNameDouble::ND_ELT_MOLES; + it_ec = ec_vector.begin(); + it_f = f_vector.begin(); + for (; it_ec != ec_vector.end(); ) { + double extensive = *it_f; + double intensive = extensive/sum_extensive; + this->moles += it_ec->moles*extensive; + this->la += 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); + it_ec++; + it_f++; + } +} cxxExchComp::~cxxExchComp() { diff --git a/ExchComp.h b/ExchComp.h index dc2f7abb..64972d6d 100644 --- a/ExchComp.h +++ b/ExchComp.h @@ -18,6 +18,7 @@ class cxxExchComp public: cxxExchComp(); cxxExchComp(struct exch_comp *); + cxxExchComp(std::vector &ec_vector, std::vector &f_vector); ~cxxExchComp(); diff --git a/Exchange.h b/Exchange.h index 76306f85..c4b3e34e 100644 --- a/Exchange.h +++ b/Exchange.h @@ -35,10 +35,15 @@ public: bool get_related_rate(void); + std::list &get_exchComps(void) { + return(this->exchComps); + } + #ifdef USE_MPI - void cxxExchange::mpi_pack(std::vector& ints, std::vector& doubles); - void cxxExchange::mpi_unpack(int *ints, int *ii, double *doubles, int *dd); + void mpi_pack(std::vector& ints, std::vector& doubles); + void mpi_unpack(int *ints, int *ii, double *doubles, int *dd); #endif + protected: std::list exchComps; bool pitzer_exchange_gammas; diff --git a/NameDouble.h b/NameDouble.h index ab6c5b07..cac0a527 100644 --- a/NameDouble.h +++ b/NameDouble.h @@ -46,6 +46,10 @@ public: void add(const cxxNameDouble &old, double factor); + void insert(char *str, double d) { + (*this)[str] = d; + } + void mpi_pack(std::vector& ints, std::vector& doubles); void mpi_pack(int *ints, int *ii, double *doubles, int *dd); diff --git a/Solution.cxx b/Solution.cxx index c47bfebc..3f602d52 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -1003,6 +1003,50 @@ void cxxSolution::set_master_activity(char *string, double d) { this->master_activity[string] = d; } +#ifdef SKIP +#include "../hst.h" +/* ---------------------------------------------------------------------- */ +void cxxSolution::xsolution_save(int n) +/* ---------------------------------------------------------------------- */ +{ +/* + * Save solution composition into Solution class + * + * input: n is pointer number in solution + */ + this->set_description(" "); + this->n_user = n; + this->n_user_end = n; + this->tc = tc_x; + this->ph = ph_x; + this->pe = solution_pe_x; + this->mu = mu_x; + this->ah2o = ah2o_x; + this->total_h = total_h_x; + this->total_o = total_o_x; + this->cb = cb_x; + this->mass_water = mass_water_aq_x; + this->total_alkalinity = total_alkalinity; +/* + * Copy totals data + */ + for (int j = 2; j < count_total; j++) { + this->totals.insert(buffer[j].master->elt->name, buffer[j].master->total_primary); + } + + for (int j = 0; j < count_activity_list; j++) { + this->master_activity.insert(activity_list[j].master->elt->name, activity_list[j].master->s->la); + } + if (pitzer_model == TRUE) { + for (int j= 0; j < count_s; j++) { + if (s[j]->lg != 0.0) { + this->species_gamma.insert(s[j]->name, s[j]->lg); + } + } + } +} +#endif + #include "ISolution.h" #include "Exchange.h" diff --git a/StorageBin.cxx b/StorageBin.cxx index fe450649..71c0016d 100644 --- a/StorageBin.cxx +++ b/StorageBin.cxx @@ -321,7 +321,6 @@ void cxxStorageBin::add(struct system * system_ptr) } } -//#include // std::cout std::cerr void cxxStorageBin::cxxStorageBin2phreeqc(int n) // // copy data fromphreeqc storage to storage bin @@ -818,3 +817,85 @@ void cxxStorageBin::mpi_recv(int task_number) } #endif + +cxxExchange *cxxStorageBin::mix_cxxExchange(cxxMix &mixmap) + +{ +/* + * mixes exchangers based on cxxMix structure, returns new exchanger + * return exchanger must be freed by calling method + */ + cxxExchange *new_exch_ptr, *old_exch_ptr; +/* + * Zero out global solution data + */ + new_exch_ptr = new cxxExchange(); + + std::map::const_iterator it_mix; + std::map *mixcomps = mixmap.comps(); +/* + * Make list of ExchComps + */ + std::vector ec_vector; + std::vector f_vector; + // + // make list of all exchange components and their mix fractions + // + for (it_mix = mixcomps->begin(); it_mix != mixcomps->end(); it_mix++) { + old_exch_ptr = &((this->Exchangers.find(it_mix->first))->second); + if (old_exch_ptr == NULL) { + sprintf(error_string, "Exchange %d not found in mix_cxxExchange.", it_mix->first); + error_msg(error_string, CONTINUE); + input_error++; + return(NULL); + } + // Add exchange components to vector ec_vector + std::list::const_iterator it_ec; + std::list &eclist = old_exch_ptr->get_exchComps(); + for (it_ec = eclist.begin(); it_ec != eclist.end(); it_ec++) { + f_vector.push_back(it_mix->second); + //cxxExchComp ec = *it_ec; + //ec_vector.push_back(ec); + ec_vector.push_back(*it_ec); + } + } + // + // Process list to make mixture + // + char *current_formula = ec_vector.begin()->get_formula(); + while (current_formula != NULL) { + + std::vector ec_subvector; + std::vector f_subvector; + std::vector::iterator it_ec = ec_vector.begin(); + std::vector::iterator it_f = f_vector.begin(); + current_formula = NULL; + for ( ; it_ec != ec_vector.end(); it_ec++) { + if (*it_f != 0) { + if (current_formula == NULL) current_formula = it_ec->get_formula(); + if (it_ec->get_formula() == current_formula) { + ec_subvector.push_back(*it_ec); + f_subvector.push_back(*it_f); + *it_f = 0; + //ec_vector.erase(it_ec); + //f_vector.erase(it_f); + } + } + it_f++; + } + // + // mix ec_subvector to make + // one exchange component + // + if (current_formula != NULL) { + cxxExchComp new_comp(ec_subvector, f_subvector); + new_exch_ptr->get_exchComps().push_back(new_comp); + } + } + /* + std::ostringstream oss; + new_exch_ptr->dump_raw(oss, 0); + std::cerr << oss.str(); + */ + return(new_exch_ptr); +} diff --git a/StorageBin.h b/StorageBin.h index c65e8c2f..3b19c292 100644 --- a/StorageBin.h +++ b/StorageBin.h @@ -138,6 +138,7 @@ public: struct system *cxxStorageBin2system(int i); cxxSolution *mix_cxxSolutions(cxxMix &mixmap); + cxxExchange *mix_cxxExchange(cxxMix &mixmap); #ifdef USE_MPI void mpi_send(int n, int task_number);