From 566bfbce3bc5a8616d57f6e26a927aa94f44183c Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sun, 21 Feb 2021 21:13:05 -0700 Subject: [PATCH 1/4] Added GetSpeciesLog10Molalities. Tested OpenMP with VS. Tested MPI with MinGW. Fortran, C, and C++ seem to work. --- phreeqcpp/Solution.cxx | 80 ++++++++++++++++++++++++++++++++++++++++-- phreeqcpp/Solution.h | 2 ++ phreeqcpp/mainsubs.cpp | 9 +++++ 3 files changed, 89 insertions(+), 2 deletions(-) diff --git a/phreeqcpp/Solution.cxx b/phreeqcpp/Solution.cxx index 28e11fdf..4b454ab1 100644 --- a/phreeqcpp/Solution.cxx +++ b/phreeqcpp/Solution.cxx @@ -83,6 +83,7 @@ cxxSolution::operator =(const cxxSolution &rhs) this->isotopes = rhs.isotopes; this->species_map = rhs.species_map; this->log_gamma_map = rhs.log_gamma_map; + this->log_molalities_map = rhs.log_molalities_map; if (this->initial_data) delete initial_data; if (rhs.initial_data != NULL) @@ -327,6 +328,19 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con s_oss << it->first << " " << it->second << "\n"; } } + + // log_molalities_map + if (log_molalities_map.size() > 0) + { + s_oss << indent1; + s_oss << "-log_molalities_map" << "\n"; + std::map::const_iterator it = this->log_molalities_map.begin(); + for (; it != log_molalities_map.end(); it++) + { + s_oss << indent2; + s_oss << it->first << " " << it->second << "\n"; + } + } return; } @@ -1013,7 +1027,32 @@ cxxSolution::read_raw(CParser & parser, bool check) } opt_save = CParser::OPT_DEFAULT; break; - + case 27: // log_molalities_map + { + int s_num; + if (parser.peek_token() != CParser::TT_EMPTY) + { + if (!(parser.get_iss() >> s_num)) + { + parser.incr_input_error(); + parser.error_msg("Expected integer for species number.", + PHRQ_io::OT_CONTINUE); + } + else + { + double d; + if (!(parser.get_iss() >> d)) + { + parser.incr_input_error(); + parser.error_msg("Expected double for species molality.", + PHRQ_io::OT_CONTINUE); + } + this->log_molalities_map[s_num] = d; + } + } + opt_save = 27; + } + break; } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -1415,6 +1454,19 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive) this->log_gamma_map[git->first] = git->second; } } + // Add molalities + std::map::const_iterator mit = addee.log_molalities_map.begin(); + for (; mit != addee.log_molalities_map.end(); mit++) + { + if (this->log_molalities_map.find(mit->first) != this->log_molalities_map.end()) + { + this->log_molalities_map[mit->first] = this->log_molalities_map[mit->first] * f1 + mit->second * f2; + } + else + { + this->log_molalities_map[mit->first] = mit->second; + } + } } } @@ -1618,6 +1670,18 @@ cxxSolution::Serialize(Dictionary & dictionary, std::vector < int >&ints, doubles.push_back(it->second); } } + /* + * log_molalities_map + */ + ints.push_back((int)log_molalities_map.size()); + { + std::map < int, double >::iterator it; + for (it = log_molalities_map.begin(); it != log_molalities_map.end(); it++) + { + ints.push_back(it->first); + doubles.push_back(it->second); + } + } } /* ---------------------------------------------------------------------- */ @@ -1692,6 +1756,17 @@ cxxSolution::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std: log_gamma_map[ints[ii++]] = doubles[dd++]; } } + /* + * log_molalities_map + */ + { + log_molalities_map.clear(); + int n = ints[ii++]; + for (int i = 0; i < n; i++) + { + log_molalities_map[ints[ii++]] = doubles[dd++]; + } + } } @@ -1722,6 +1797,7 @@ const std::vector< std::string >::value_type temp_vopts[] = { std::vector< std::string >::value_type("soln_vol"), // 23 std::vector< std::string >::value_type("species_map"), // 24 std::vector< std::string >::value_type("log_gamma_map"), // 25 - std::vector< std::string >::value_type("potential") // 26 + std::vector< std::string >::value_type("potential"), // 26 + std::vector< std::string >::value_type("log_molalities_map") // 27 }; const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); \ No newline at end of file diff --git a/phreeqcpp/Solution.h b/phreeqcpp/Solution.h index 1ab18654..cecabe78 100644 --- a/phreeqcpp/Solution.h +++ b/phreeqcpp/Solution.h @@ -66,6 +66,7 @@ class cxxSolution:public cxxNumKeyword cxxNameDouble & Get_species_gamma(void) {return this->species_gamma;} std::map & Get_species_map(void) {return this->species_map;} std::map & Get_log_gamma_map(void) {return this->log_gamma_map;} + std::map& Get_log_molalities_map(void) { return this->log_molalities_map; } std::map < std::string, cxxSolutionIsotope > & Get_isotopes(void) {return this->isotopes;} const std::map < std::string, cxxSolutionIsotope > & Get_isotopes(void)const {return this->isotopes;} void Set_isotopes(const std::map < std::string, cxxSolutionIsotope > &iso ) {this->isotopes = iso;} @@ -141,6 +142,7 @@ class cxxSolution:public cxxNumKeyword const static std::vector < std::string > vopts; std::map species_map; std::map log_gamma_map; + std::map log_molalities_map; }; #endif // !defined(SOLUTION_H_INCLUDED) diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 27e34981..efe65482 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -1671,6 +1671,15 @@ xsolution_save(int n_user) temp_solution.Get_log_gamma_map()[s_x[i]->number] = s_x[i]->lg; } } + // saves molalities + temp_solution.Get_log_molalities_map().clear(); + for (int i = 0; i < this->count_s_x; i++) + { + if (s_x[i]->type <= H2O) + { + temp_solution.Get_log_molalities_map()[s_x[i]->number] = s_x[i]->lm; + } + } } /* * Save solution From 0a63328909ececf7f6b15d7136d9460883f0768f Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Tue, 23 Feb 2021 10:56:26 -0700 Subject: [PATCH 2/4] C++ is working with OpenMP and MPI for Get/SetGasPhaseMoles. Need to add c and F90. --- phreeqcpp/GasPhase.cxx | 48 ++++++++++++++++++++++++++++++++++++++++++ phreeqcpp/GasPhase.h | 4 +++- 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/phreeqcpp/GasPhase.cxx b/phreeqcpp/GasPhase.cxx index 7a538593..68c82adf 100644 --- a/phreeqcpp/GasPhase.cxx +++ b/phreeqcpp/GasPhase.cxx @@ -649,6 +649,54 @@ LDBLE cxxGasPhase::Calc_total_moles(void)const } return tot; } +void cxxGasPhase::Delete_component(const std::string comp_name) +{ + for (size_t i = 0; i < this->gas_comps.size(); i++) + { + if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0) + { + this->gas_comps.erase(this->gas_comps.begin() + i); // To delete the ith element + break; + } + } +} +void cxxGasPhase::Set_component_moles(const std::string comp_name, const double moles) +{ + size_t i; + if (moles < 0.0) + { + this->Delete_component(comp_name); + } + else + { + cxxGasComp* ptr = this->Find_comp(comp_name.c_str()); + if (ptr != NULL) + { + ptr->Set_moles(moles); + } + else + { + cxxGasComp temp_comp; + temp_comp.Set_phase_name(comp_name); + temp_comp.Set_moles(moles); + this->gas_comps.push_back(temp_comp); + } + } +} +double cxxGasPhase::Get_component_moles(const std::string comp_name) +{ + double moles = -1.0; + for (size_t i = 0; i < this->gas_comps.size(); i++) + { + if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0) + { + moles = this->gas_comps[i].Get_moles(); + break; + } + } + return moles; +} + cxxGasComp * cxxGasPhase::Find_comp(const char * comp_name) { diff --git a/phreeqcpp/GasPhase.h b/phreeqcpp/GasPhase.h index 5741f10c..9da519dd 100644 --- a/phreeqcpp/GasPhase.h +++ b/phreeqcpp/GasPhase.h @@ -70,7 +70,9 @@ class cxxGasPhase:public cxxNumKeyword cxxGasComp *Find_comp(const char * comp_name); void Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles); void Deserialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles, int &ii, int &dd); - + void Delete_component(const std::string comp_name); + void Set_component_moles(const std::string comp_name, const double moles); + double Get_component_moles(const std::string comp_name); protected: void add(const cxxGasPhase & addee, LDBLE extensive); From e1554d9987f2b234df485aabcfa44de18641616c Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Thu, 25 Feb 2021 18:04:29 -0700 Subject: [PATCH 3/4] more checking in. Should be down to tweaks for SetGasPhaseMoles. --- phreeqcpp/GasPhase.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/phreeqcpp/GasPhase.cxx b/phreeqcpp/GasPhase.cxx index 68c82adf..1d694f75 100644 --- a/phreeqcpp/GasPhase.cxx +++ b/phreeqcpp/GasPhase.cxx @@ -662,7 +662,6 @@ void cxxGasPhase::Delete_component(const std::string comp_name) } void cxxGasPhase::Set_component_moles(const std::string comp_name, const double moles) { - size_t i; if (moles < 0.0) { this->Delete_component(comp_name); From 18158db399c6d6ad6d77b24912559deee0e78c85 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 27 Feb 2021 22:26:48 -0700 Subject: [PATCH 4/4] Last of changes for GetGasPhasePressures and GetGasPhasePhi, openmp and mpi. MPI fortrans not tested. --- phreeqcpp/GasComp.cxx | 71 ++++++++++++++++++++++++++++++++++++------ phreeqcpp/GasComp.h | 9 ++++++ phreeqcpp/GasPhase.cxx | 39 +++++++++++++++++++++++ phreeqcpp/GasPhase.h | 3 ++ phreeqcpp/mainsubs.cpp | 46 ++++++++++++++++++--------- phreeqcpp/tidy.cpp | 36 ++++++++++++++++----- 6 files changed, 172 insertions(+), 32 deletions(-) diff --git a/phreeqcpp/GasComp.cxx b/phreeqcpp/GasComp.cxx index 21a13c71..59f3d796 100644 --- a/phreeqcpp/GasComp.cxx +++ b/phreeqcpp/GasComp.cxx @@ -30,6 +30,9 @@ cxxGasComp::cxxGasComp(PHRQ_io *io) p_read = 0.0; moles = 0.0; initial_moles = 0.0; + p = 0.0; + phi = 0.0; + f = 0.0; } cxxGasComp::~cxxGasComp() @@ -56,6 +59,9 @@ cxxGasComp::dump_raw(std::ostream & s_oss, unsigned int indent) const s_oss << indent0 << "# GasComp workspace variables #\n"; s_oss << indent0 << "-initial_moles " << this->initial_moles << "\n"; + s_oss << indent0 << "-p " << this->p << "\n"; + s_oss << indent0 << "-phi " << this->phi << "\n"; + s_oss << indent0 << "-f " << this->f << "\n"; } bool @@ -125,6 +131,37 @@ cxxGasComp::read_raw(CParser & parser, bool check) PHRQ_io::OT_CONTINUE); } break; + + case 5: // p + if (!(parser.get_iss() >> this->p)) + { + this->p = 0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for pressure.", + PHRQ_io::OT_CONTINUE); + } + break; + + case 6: // phi + if (!(parser.get_iss() >> this->phi)) + { + this->phi = 0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for phi.", + PHRQ_io::OT_CONTINUE); + } + break; + + case 7: // f + if (!(parser.get_iss() >> this->f)) + { + this->f = 0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for f.", + PHRQ_io::OT_CONTINUE); + } + break; + } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -150,9 +187,9 @@ cxxGasComp::add(const cxxGasComp & addee, LDBLE extensive) if (addee.phase_name.size() == 0) return; - /* - ext1 = this->moles; - ext2 = addee.moles * extensive; + double f1, f2; + double ext1 = this->moles; + double ext2 = addee.moles * extensive; if (ext1 + ext2 != 0) { f1 = ext1 / (ext1 + ext2); @@ -163,13 +200,15 @@ cxxGasComp::add(const cxxGasComp & addee, LDBLE extensive) f1 = 0.5; f2 = 0.5; } - */ assert(this->phase_name == addee.phase_name); - this->p_read += addee.p_read * extensive; + this->p_read = this->p_read*f1 + addee.p_read * f2; this->moles += addee.moles * extensive; this->initial_moles += addee.initial_moles * extensive; + this->p = this->p * f1 + addee.p * f2; + this->phi = this->phi * f1 + addee.phi * f2; + this->f = this->f * f1 + addee.f * f2; } void @@ -178,6 +217,9 @@ cxxGasComp::multiply(LDBLE extensive) this->p_read *= extensive; this->moles *= extensive; this->initial_moles *= extensive; + this->p *= 1.0; + this->phi *= 1.0; + this->f *= 1.0; } void cxxGasComp::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles) @@ -186,6 +228,9 @@ cxxGasComp::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::ve doubles.push_back(this->moles); doubles.push_back(this->p_read); doubles.push_back(this->initial_moles); + doubles.push_back(this->p); + doubles.push_back(this->phi); + doubles.push_back(this->f); } void @@ -196,13 +241,19 @@ cxxGasComp::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->moles = doubles[dd++]; this->p_read = doubles[dd++]; this->initial_moles = doubles[dd++]; + this->p = doubles[dd++]; + this->phi = doubles[dd++]; + this->f = doubles[dd++]; } const std::vector< std::string >::value_type temp_vopts[] = { - std::vector< std::string >::value_type("phase_name"), // 0 - std::vector< std::string >::value_type("name"), // 1 - std::vector< std::string >::value_type("p_read"), // 2 - std::vector< std::string >::value_type("moles"), // 3 - std::vector< std::string >::value_type("initial_moles") // 4 + std::vector< std::string >::value_type("phase_name"), // 0 + std::vector< std::string >::value_type("name"), // 1 + std::vector< std::string >::value_type("p_read"), // 2 + std::vector< std::string >::value_type("moles"), // 3 + std::vector< std::string >::value_type("initial_moles"), // 4 + std::vector< std::string >::value_type("p"), // 5 + std::vector< std::string >::value_type("phi"), // 6 + std::vector< std::string >::value_type("f") // 7 }; const std::vector< std::string > cxxGasComp::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/phreeqcpp/GasComp.h b/phreeqcpp/GasComp.h index d4661a6b..ce44e100 100644 --- a/phreeqcpp/GasComp.h +++ b/phreeqcpp/GasComp.h @@ -28,6 +28,12 @@ class cxxGasComp: public PHRQ_base void Set_p_read(LDBLE t) {this->p_read = t;} LDBLE Get_moles() const {return this->moles;} void Set_moles(LDBLE t) {this->moles = t;} + LDBLE Get_p() const { return this->p; } + void Set_p(LDBLE t) { this->p = t; } + LDBLE Get_phi() const { return this->phi; } + void Set_phi(LDBLE t) { this->phi = t; } + LDBLE Get_f() const { return this->f; } + void Set_f(LDBLE t) { this->f = t; } LDBLE Get_initial_moles() const {return this->initial_moles;} void Set_initial_moles(LDBLE t) {this->initial_moles = t;} @@ -44,6 +50,9 @@ class cxxGasComp: public PHRQ_base LDBLE p_read; // internal workspace LDBLE initial_moles; + LDBLE p; + LDBLE phi; + LDBLE f; const static std::vector < std::string > vopts; }; diff --git a/phreeqcpp/GasPhase.cxx b/phreeqcpp/GasPhase.cxx index 1d694f75..ad992501 100644 --- a/phreeqcpp/GasPhase.cxx +++ b/phreeqcpp/GasPhase.cxx @@ -695,6 +695,45 @@ double cxxGasPhase::Get_component_moles(const std::string comp_name) } return moles; } +double cxxGasPhase::Get_component_p(const std::string comp_name) +{ + double p = -1.0; + for (size_t i = 0; i < this->gas_comps.size(); i++) + { + if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0) + { + p = this->gas_comps[i].Get_p(); + break; + } + } + return p; +} +double cxxGasPhase::Get_component_phi(const std::string comp_name) +{ + double phi = -1.0; + for (size_t i = 0; i < this->gas_comps.size(); i++) + { + if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0) + { + phi = this->gas_comps[i].Get_phi(); + break; + } + } + return phi; +} +double cxxGasPhase::Get_component_f(const std::string comp_name) +{ + double f = -1.0; + for (size_t i = 0; i < this->gas_comps.size(); i++) + { + if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0) + { + f = this->gas_comps[i].Get_f(); + break; + } + } + return f; +} cxxGasComp * cxxGasPhase::Find_comp(const char * comp_name) diff --git a/phreeqcpp/GasPhase.h b/phreeqcpp/GasPhase.h index 9da519dd..1798feb1 100644 --- a/phreeqcpp/GasPhase.h +++ b/phreeqcpp/GasPhase.h @@ -73,6 +73,9 @@ class cxxGasPhase:public cxxNumKeyword void Delete_component(const std::string comp_name); void Set_component_moles(const std::string comp_name, const double moles); double Get_component_moles(const std::string comp_name); + double Get_component_p(const std::string comp_name); + double Get_component_phi(const std::string comp_name); + double Get_component_f(const std::string comp_name); protected: void add(const cxxGasPhase & addee, LDBLE extensive); diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 27e34981..e70726b0 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -1352,19 +1352,19 @@ int Phreeqc:: xgas_save(int n_user) /* ---------------------------------------------------------------------- */ { -/* - * Save gas composition into structure gas_phase with user - * number n_user. - */ + /* + * Save gas composition into structure gas_phase with user + * number n_user. + */ char token[MAX_LENGTH]; if (use.Get_gas_phase_ptr() == NULL) return (OK); - cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr(); + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); cxxGasPhase temp_gas_phase(*gas_phase_ptr); -/* - * Store in gas_phase - */ + /* + * Store in gas_phase + */ temp_gas_phase.Set_n_user(n_user); temp_gas_phase.Set_n_user_end(n_user); sprintf(token, "Gas phase after simulation %d.", simulation); @@ -1372,22 +1372,38 @@ xgas_save(int n_user) temp_gas_phase.Set_new_def(false); temp_gas_phase.Set_solution_equilibria(false); temp_gas_phase.Set_n_solution(-99); -/* - * Update amounts - */ - for (size_t i = 0 ; i < temp_gas_phase.Get_gas_comps().size(); i++) + /* + * Update amounts + */ + bool PR = false; + if (gas_phase_ptr->Get_v_m() >= 0.01) PR = true; + for (size_t i = 0; i < temp_gas_phase.Get_gas_comps().size(); i++) { - cxxGasComp * gc_ptr = &(temp_gas_phase.Get_gas_comps()[i]); + cxxGasComp* gc_ptr = &(temp_gas_phase.Get_gas_comps()[i]); int k; - struct phase *phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE); + struct phase* phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE); assert(phase_ptr); - gc_ptr->Set_moles(phase_ptr->moles_x); + if (PR) + { + gc_ptr->Set_moles(phase_ptr->moles_x); + gc_ptr->Set_p(phase_ptr->p_soln_x); + gc_ptr->Set_phi(phase_ptr->pr_phi); + gc_ptr->Set_f(phase_ptr->p_soln_x * phase_ptr->pr_phi); + } + else + { + gc_ptr->Set_moles(phase_ptr->moles_x); + gc_ptr->Set_p(phase_ptr->p_soln_x); + gc_ptr->Set_phi(1.0); + gc_ptr->Set_f(phase_ptr->p_soln_x); + } } Rxn_gas_phase_map[n_user] = temp_gas_phase; use.Set_gas_phase_ptr(NULL); return (OK); } + /* ---------------------------------------------------------------------- */ int Phreeqc:: xss_assemblage_save(int n_user) diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp index 95d12e21..152af18b 100644 --- a/phreeqcpp/tidy.cpp +++ b/phreeqcpp/tidy.cpp @@ -971,9 +971,14 @@ tidy_gas_phase(void) { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); if (!PR) - gas_phase_ptr->Get_gas_comps()[j].Set_moles( - gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * gas_phase_ptr->Get_volume() / - R_LITER_ATM / gas_phase_ptr->Get_temperature()); + { + double moles = gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * gas_phase_ptr->Get_volume() / + R_LITER_ATM / gas_phase_ptr->Get_temperature(); + gas_phase_ptr->Get_gas_comps()[j].Set_moles(moles); + gas_phase_ptr->Get_gas_comps()[j].Set_p(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()); + gas_phase_ptr->Get_gas_comps()[j].Set_phi(1.0); + gas_phase_ptr->Get_gas_comps()[j].Set_f(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()); + } } else { @@ -999,10 +1004,15 @@ tidy_gas_phase(void) { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); if (!PR) - gas_phase_ptr->Get_gas_comps()[j].Set_moles ( - gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * - gas_phase_ptr->Get_volume() / R_LITER_ATM / - gas_phase_ptr->Get_temperature()); + { + double moles = gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * + gas_phase_ptr->Get_volume() / R_LITER_ATM / + gas_phase_ptr->Get_temperature(); + gas_phase_ptr->Get_gas_comps()[j].Set_moles(moles); + gas_phase_ptr->Get_gas_comps()[j].Set_p(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()); + gas_phase_ptr->Get_gas_comps()[j].Set_phi(1.0); + gas_phase_ptr->Get_gas_comps()[j].Set_f(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()); + } } else { @@ -1027,7 +1037,13 @@ tidy_gas_phase(void) int k; struct phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j_PR].Get_phase_name().c_str(), &k, FALSE); if (gc[j_PR].Get_p_read() == 0) + { + gc[j_PR].Set_moles(0.0); + gc[j_PR].Set_p(0.0); + gc[j_PR].Set_phi(1.0); + gc[j_PR].Set_f(0.0); continue; + } if (phase_ptr) { phase_ptr->moles_x = gc[j_PR].Get_p_read() / P; @@ -1047,11 +1063,17 @@ tidy_gas_phase(void) if (gc[j_PR].Get_p_read() == 0) { gc[j_PR].Set_moles(0.0); + gc[j_PR].Set_p(0.0); + gc[j_PR].Set_phi(1.0); + gc[j_PR].Set_f(0.0); } else { if (phase_ptr) { gc[j_PR].Set_moles(phase_ptr->moles_x * gas_phase_ptr->Get_volume() / V_m); + gc[j_PR].Set_p(gc[j_PR].Get_p_read()); + gc[j_PR].Set_phi(phase_ptr->pr_phi); + gc[j_PR].Set_f(gc[j_PR].Get_p_read()* phase_ptr->pr_phi); gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gc[j_PR].Get_moles()); } }