diff --git a/phreeqcpp/ChartObject.cpp b/phreeqcpp/ChartObject.cpp index 9d6513ce..018e3b25 100644 --- a/phreeqcpp/ChartObject.cpp +++ b/phreeqcpp/ChartObject.cpp @@ -11,7 +11,8 @@ #include "ChartObject.h" #include "Parser.h" #include -#include +//#include +#include #include #include "phqalloc.h" diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index 99b4661d..b76897ad 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -7215,7 +7215,7 @@ _NilCheck(void) return _Escape(-3); } -#ifdef SKIP +#ifdef NPP /* The following is suitable for the HP Pascal operating system. It might want to be revised when emulating another system. */ diff --git a/phreeqcpp/PBasic.h b/phreeqcpp/PBasic.h index 95d77db2..73321cc8 100644 --- a/phreeqcpp/PBasic.h +++ b/phreeqcpp/PBasic.h @@ -7,7 +7,8 @@ #include #include #include -#include +//#include +#include #include #include "phrqtype.h" #include "PHRQ_base.h" diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index a0e38bc0..9881c015 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -283,12 +283,12 @@ public: int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1); int calc_all_donnan(void); int calc_init_donnan(void); + LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector &zcorr); LDBLE g_function(LDBLE x_value); LDBLE midpnt(LDBLE x1, LDBLE x2, int n); void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv, LDBLE* dy); LDBLE qromb_midpnt(cxxSurfaceCharge* charge_ptr, LDBLE x1, LDBLE x2); - LDBLE calc_psi_avg(cxxSurfaceCharge* charge_ptr, LDBLE surf_chrg_eq); // inverse.cpp ------------------------------- int inverse_models(void); @@ -1851,7 +1851,7 @@ isfinite handling #elif defined(HAVE_FINITE) # define PHR_ISFINITE(x) finite(x) #elif defined(HAVE_ISNAN) -# define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!isnan(x)) && ((x) != (2.0 * (x)))) ) +# define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!std::isnan(x)) && ((x) != (2.0 * (x)))) ) #else # define PHR_ISFINITE(x) ( ((x) == 0.0) || (((x) == (x)) && ((x) != (2.0 * (x)))) ) #endif diff --git a/phreeqcpp/Solution.cxx b/phreeqcpp/Solution.cxx index cb2d2941..bc1bd80e 100644 --- a/phreeqcpp/Solution.cxx +++ b/phreeqcpp/Solution.cxx @@ -47,6 +47,7 @@ cxxSolution::cxxSolution(PHRQ_io * io) this->total_o = 55.55; this->cb = 0.0; this->density = 1.0; + this->viscosity = 1.0; this->mass_water = 1.0; this->soln_vol = 1.0; this->total_alkalinity = 0.0; @@ -80,6 +81,7 @@ cxxSolution::operator =(const cxxSolution &rhs) this->total_h = rhs.total_h; this->total_o = rhs.total_o; this->density = rhs.density; + this->viscosity = rhs.viscosity; this->cb = rhs.cb; this->mass_water = rhs.mass_water; this->soln_vol = rhs.soln_vol; @@ -269,6 +271,10 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con s_oss << indent1; s_oss << "-density " << this->density << "\n"; + // new identifier + s_oss << indent1; + s_oss << "-viscosity " << this->viscosity << "\n"; + // soln_total conc structures s_oss << indent1; s_oss << "-totals" << "\n"; @@ -1070,6 +1076,16 @@ cxxSolution::read_raw(CParser & parser, bool check) opt_save = 27; } break; + case 28: // viscosity + if (!(parser.get_iss() >> this->viscosity)) + { + this->viscosity = 1.0; + parser.incr_input_error(); + parser.error_msg("Expected numeric value for viscosity.", + PHRQ_io::OT_CONTINUE); + } + opt_save = CParser::OPT_DEFAULT; + break; } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -1365,6 +1381,7 @@ cxxSolution::zero() this->total_o = 0.0; this->cb = 0.0; this->density = 1.0; + this->viscosity = 1.0; this->mass_water = 0.0; this->soln_vol = 0.0; this->total_alkalinity = 0.0; @@ -1397,6 +1414,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive) this->total_o += addee.total_o * extensive; this->cb += addee.cb * extensive; this->density = f1 * this->density + f2 * addee.density; + this->viscosity = f1 * this->viscosity + f2 * addee.viscosity; this->patm = f1 * this->patm + f2 * addee.patm; // this->potV = f1 * this->potV + f2 * addee.potV; // appt this->mass_water += addee.mass_water * extensive; @@ -1574,6 +1592,7 @@ cxxSolution::Serialize(Dictionary & dictionary, std::vector < int >&ints, doubles.push_back(this->cb); doubles.push_back(this->mass_water); doubles.push_back(this->density); + doubles.push_back(this->viscosity); doubles.push_back(this->soln_vol); doubles.push_back(this->total_alkalinity); /* @@ -1660,6 +1679,7 @@ cxxSolution::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std: this->cb = doubles[dd++]; this->mass_water = doubles[dd++]; this->density = doubles[dd++]; + this->viscosity = doubles[dd++]; this->soln_vol = doubles[dd++]; this->total_alkalinity = doubles[dd++]; /* @@ -1752,6 +1772,7 @@ const std::vector< std::string >::value_type temp_vopts[] = { 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("log_molalities_map") // 27 + std::vector< std::string >::value_type("log_molalities_map"), // 27 + std::vector< std::string >::value_type("viscosity") // 28 }; const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/phreeqcpp/Solution.h b/phreeqcpp/Solution.h index cecabe78..eef728d4 100644 --- a/phreeqcpp/Solution.h +++ b/phreeqcpp/Solution.h @@ -49,6 +49,8 @@ class cxxSolution:public cxxNumKeyword void Set_cb(LDBLE l_cb) {this->cb = l_cb;} LDBLE Get_density() const {return this->density;} void Set_density(LDBLE l_density) {this->density = l_density;} + LDBLE Get_viscosity() const { return this->viscosity; } + void Set_viscosity(LDBLE l_viscos) { this->viscosity = l_viscos; } LDBLE Get_mass_water() const {return this->mass_water;} void Set_mass_water(LDBLE l_mass_water) {this->mass_water = l_mass_water;} LDBLE Get_total_alkalinity() const {return this->total_alkalinity;} @@ -131,6 +133,7 @@ class cxxSolution:public cxxNumKeyword LDBLE cb; LDBLE mass_water; LDBLE density; + LDBLE viscosity; LDBLE soln_vol; LDBLE total_alkalinity; cxxNameDouble totals; diff --git a/phreeqcpp/SolutionIsotope.cxx b/phreeqcpp/SolutionIsotope.cxx index 78d832a5..4059b0b9 100644 --- a/phreeqcpp/SolutionIsotope.cxx +++ b/phreeqcpp/SolutionIsotope.cxx @@ -70,7 +70,8 @@ cxxSolutionIsotope::dump_xml(std::ostream & s_oss, unsigned int indent) const #ifdef NPP if (!isnan(this->ratio_uncertainty)) #else - if (this->ratio_uncertainty != NAN) + //if (this->ratio_uncertainty != NAN) + if (!std::isnan(this->ratio_uncertainty)) #endif { s_oss << indent1; diff --git a/phreeqcpp/Surface.cxx b/phreeqcpp/Surface.cxx index 40c8cc0f..ed55bb2d 100644 --- a/phreeqcpp/Surface.cxx +++ b/phreeqcpp/Surface.cxx @@ -36,6 +36,7 @@ cxxSurface::cxxSurface(PHRQ_io *io) dl_type = NO_DL; sites_units = SITES_ABSOLUTE; only_counter_ions = false; + correct_GC = false; thickness = 1e-8; debye_lengths = 0.0; DDL_viscosity = 1.0; @@ -55,6 +56,7 @@ cxxNumKeyword(io) dl_type = NO_DL; sites_units = SITES_ABSOLUTE; only_counter_ions = false; + correct_GC = false; thickness = 1e-8; debye_lengths = 0.0; DDL_viscosity = 1.0; @@ -128,6 +130,8 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons s_oss << indent1; s_oss << "-only_counter_ions " << this->only_counter_ions << "\n"; s_oss << indent1; + s_oss << "-correct_GC " << this->correct_GC << "\n"; + s_oss << indent1; s_oss << "-thickness " << this->thickness << "\n"; s_oss << indent1; s_oss << "-debye_lengths " << this->debye_lengths << "\n"; @@ -189,6 +193,7 @@ cxxSurface::read_raw(CParser & parser, bool check) this->Set_tidied(true); bool only_counter_ions_defined(false); + //bool correct_GC_defined(false); bool thickness_defined(false); bool type_defined(false); bool dl_type_defined(false); @@ -468,6 +473,17 @@ cxxSurface::read_raw(CParser & parser, bool check) PHRQ_io::OT_CONTINUE); } break; + case 19: // correct_GC + if (!(parser.get_iss() >> this->correct_GC)) + { + this->correct_GC = false; + parser.incr_input_error(); + parser. + error_msg("Expected boolean value for correct_GC.", + PHRQ_io::OT_CONTINUE); + } + //correct_GC_defined = true; + break; } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -482,6 +498,13 @@ cxxSurface::read_raw(CParser & parser, bool check) error_msg("Only_counter_ions not defined for SURFACE_RAW input.", PHRQ_io::OT_CONTINUE); } + //if (correct_GC_defined == false) + //{ + // parser.incr_input_error(); + // parser. + // error_msg("correct_GC not defined for SURFACE_RAW input.", + // PHRQ_io::OT_CONTINUE); + //} if (thickness_defined == false) { parser.incr_input_error(); @@ -559,6 +582,7 @@ cxxSurface::add(const cxxSurface & addee_in, LDBLE extensive) if (this->surface_comps.size() == 0) { this->only_counter_ions = addee.only_counter_ions; + this->correct_GC = addee.correct_GC; this->dl_type = addee.dl_type; this->type = addee.type; this->sites_units = addee.sites_units; @@ -730,6 +754,7 @@ cxxSurface::Serialize(Dictionary & dictionary, std::vector < int >&ints, doubles.push_back(this->debye_lengths); doubles.push_back(this->DDL_viscosity); doubles.push_back(this->DDL_limit); + ints.push_back(this->correct_GC ? 1 : 0); ints.push_back(this->transport ? 1 : 0); this->totals.Serialize(dictionary, ints, doubles); ints.push_back(this->solution_equilibria ? 1 : 0); @@ -776,6 +801,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->debye_lengths = doubles[dd++]; this->DDL_viscosity = doubles[dd++]; this->DDL_limit = doubles[dd++]; + this->correct_GC = (ints[ii++] != 0); this->transport = (ints[ii++] != 0); this->totals.Deserialize(dictionary, ints, doubles, ii, dd); this->solution_equilibria = (ints[ii++] != 0); @@ -803,6 +829,7 @@ const std::vector< std::string >::value_type temp_vopts[] = { std::vector< std::string >::value_type("solution_equilibria"), // 15 std::vector< std::string >::value_type("n_solution"), // 16 std::vector< std::string >::value_type("totals"), // 17 - std::vector< std::string >::value_type("tidied") // 18 -}; + std::vector< std::string >::value_type("tidied"), // 18 + std::vector< std::string >::value_type("correct_gc") // 19 +}; const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/phreeqcpp/Surface.h b/phreeqcpp/Surface.h index 7f97bce3..f01c716c 100644 --- a/phreeqcpp/Surface.h +++ b/phreeqcpp/Surface.h @@ -67,6 +67,9 @@ public: void Set_DDL_viscosity(LDBLE t) {DDL_viscosity = t;} LDBLE Get_DDL_limit(void) const {return DDL_limit;} void Set_DDL_limit(LDBLE t) {DDL_limit = t;} + bool Get_correct_GC(void) const { return correct_GC; } + void Set_correct_GC(bool tf) { correct_GC = tf; } + std::vector Donnan_factors; bool Get_transport(void) const {return transport;} void Set_transport(bool tf) {transport = tf;} cxxNameDouble & Get_totals() {return this->totals;} @@ -91,6 +94,7 @@ protected: LDBLE debye_lengths; LDBLE DDL_viscosity; LDBLE DDL_limit; + bool correct_GC; bool transport; cxxNameDouble totals; bool solution_equilibria; diff --git a/phreeqcpp/basicsubs.cpp b/phreeqcpp/basicsubs.cpp index 1a6d11df..3094c18b 100644 --- a/phreeqcpp/basicsubs.cpp +++ b/phreeqcpp/basicsubs.cpp @@ -187,21 +187,60 @@ diff_c(const char *species_name) /* ---------------------------------------------------------------------- */ { class species *s_ptr; - LDBLE g; + LDBLE ka, l_z, Dw, ff, sqrt_mu; + sqrt_mu = sqrt(mu_x); s_ptr = s_search(species_name); - if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/) + //LDBLE g; + //if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/) + //{ + // g = s_ptr->dw; + // if (s_ptr->dw_t) + // g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + // g *= viscos_0_25 / viscos * tk_x / 298.15; + //} + //else + //{ + // g = 0; + //} + //return (g); + if (s_ptr == NULL) + return(0); + if ((Dw = s_ptr->dw) == 0) { - g = s_ptr->dw; - if (s_ptr->dw_t) - g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - g *= viscos_0_25 / viscos * tk_x / 298.15; + if (correct_Dw) + Dw = default_Dw; + else + return(0); + } + if ((l_z = fabs(s_ptr->z)) == 0) + { + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field } else { - g = 0; + if (s_ptr->dw_a2) + ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); + else + ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); + if (s_ptr->dw_a) + { + ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); + //if (print_viscosity && s_ptr->dw_a_visc) + // ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc); + } + else + { + ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + } + Dw *= ff; } - return (g); + + if (tk_x != 298.15 && s_ptr->dw_t) + Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + + s_ptr->dw_corr = Dw; + return (Dw * viscos_0_25 / viscos_0); } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: @@ -209,22 +248,58 @@ setdiff_c(const char *species_name, double d) /* ---------------------------------------------------------------------- */ { class species *s_ptr; - LDBLE g; + LDBLE ka, l_z, Dw, ff, sqrt_mu; + sqrt_mu = sqrt(mu_x); s_ptr = s_search(species_name); - if (s_ptr != NULL) + + //LDBLE g; + //s_ptr = s_search(species_name); + //if (s_ptr != NULL) + //{ + + // s_ptr->dw = d; + // g = s_ptr->dw; + // if (s_ptr->dw_t) + // g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + // g *= viscos_0_25 / viscos * tk_x / 298.15;; + //} + //else + //{ + // g = 0; + //} + //return (g); + if (s_ptr == NULL) + return(0); + Dw = s_ptr->dw = d; + if ((l_z = fabs(s_ptr->z)) == 0) { - s_ptr->dw = d; - g = s_ptr->dw; - if (s_ptr->dw_t) - g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); - g *= viscos_0_25 / viscos * tk_x / 298.15;; + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field } else { - g = 0; + if (s_ptr->dw_a2) + ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); + else + ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); + if (s_ptr->dw_a) + { + ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); + //if (print_viscosity && s_ptr->dw_a_visc) + // ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc); + } + else + { + ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + } + Dw *= ff; } - return (g); + + if (tk_x != 298.15 && s_ptr->dw_t) + Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15); + + s_ptr->dw_corr = Dw; + return (Dw * viscos_0_25 / viscos_0); } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: @@ -314,19 +389,18 @@ calc_SC(void) SC = 0; //LDBLE ta1, ta2, ta3, ta4; - for (i = 0; i < (int)this->s_x.size(); i++) - { - // ** for optimizing, get numbers from -analyt for H+ = H+... - //if (!strcmp(s_x[i]->name, "H+")) - //{ - // ta1 = s_x[i]->logk[2]; - // ta2 = s_x[i]->logk[3]; - // ta3 = s_x[i]->logk[4]; - // ta4 = s_x[i]->logk[5]; - // break; - //} - // - } + //for (i = 0; i < (int)this->s_x.size(); i++) + //{ + // // ** for optimizing, get numbers from -analyt for H+ = H+... + // if (!strcmp(s_x[i]->name, "H+")) + // { + // ta1 = s_x[i]->logk[2] * 1e15; + // ta2 = s_x[i]->logk[3] * 1e15; + // ta3 = s_x[i]->logk[4] * 1e15; + // ta4 = s_x[i]->logk[5] * 1e15; + // break; + // } + //} for (i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type != AQ && s_x[i]->type != HPLUS) @@ -338,43 +412,51 @@ calc_SC(void) else continue; } + if (s_x[i]->lm < min_dif_LM) + continue; if ((l_z = fabs(s_x[i]->z)) == 0) - l_z = 1; // only a 1st approximation for correct_Dw in electrical field - if (s_x[i]->dw_t) - Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity multiplier is done in SC - if (s_x[i]->dw_a2) - ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); + { + //l_z = 1; // only a 1st approximation for correct_Dw in electrical field + } else { - ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x , 0.75)); - //ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2)); - //ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2); - } - if (s_x[i]->dw_a) - { - ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); - if (print_viscosity) + if (s_x[i]->dw_a2) + ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75)); + else { - ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc); + ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75)); + //ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2)); + //ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2); } + if (s_x[i]->dw_a) + { + ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka)); + //if (print_viscosity && s_x[i]->dw_a_visc) + // ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc); + } + else + { + ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); + //ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka)); + } + Dw *= ff; } - else + if (tk_x != 298.15) { - ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka)); - //ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka)); + if (s_x[i]->dw_t) + Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); + //else + //{ + // Dw *= exp(ta1 / tk_x - ta1 / 298.15); + //} } - Dw *= ff; s_x[i]->dw_corr = Dw; - if (s_x[i]->z == 0) - continue; s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw; SC += s_x[i]->dw_t_SC; } SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0); - /* correct for temperature dependency... - SC_T = SC_298 * (Dw_T / T) * (298 / Dw_298) and - Dw_T = Dw_298 * (T / 298) * (viscos_298 / viscos_T) give: + /* correct for viscosity dependency... SC_T = SC_298 * (viscos_298 / viscos_T) */ SC *= viscos_0_25 / viscos_0; @@ -1520,7 +1602,8 @@ get_calculate_value(const char *name) #ifdef NPP if (isnan(rate_moles)) #else - if (rate_moles == NAN) + //if (rate_moles == NAN) + if(std::isnan(rate_moles)) #endif { error_string = sformatf( "Calculated value not SAVEed for %s.", @@ -1721,22 +1804,22 @@ pressure(void) /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -pr_phi(const char* phase_name) +pr_phi(const char *phase_name) /* ---------------------------------------------------------------------- */ { cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); if (gas_phase_ptr != NULL) { - int l; + int l; class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE); - if (phase_ptr == NULL) - { - error_string = sformatf("Gas %s, not found.", phase_name); - warning_msg(error_string); - return (1e-99); - } + if (phase_ptr == NULL) + { + error_string = sformatf( "Gas %s, not found.", phase_name); + warning_msg(error_string); + return (1e-99); + } for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++) - { + { const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]); int j; class phase* phase_ptr_gas = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE); @@ -1744,8 +1827,8 @@ pr_phi(const char* phase_name) { if (gas_phase_ptr->Get_pr_in()) { - return phase_ptr->pr_phi; - } + return phase_ptr->pr_phi; + } else { return gas_comp_ptr->Get_phi(); @@ -1753,7 +1836,7 @@ pr_phi(const char* phase_name) } } } - return(1.0); + return (1.0); } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: diff --git a/phreeqcpp/cl1.cpp b/phreeqcpp/cl1.cpp index b0dd1be4..4a3a8bbc 100644 --- a/phreeqcpp/cl1.cpp +++ b/phreeqcpp/cl1.cpp @@ -1,6 +1,7 @@ #include #include -#include +//#include +#include #include #include "Phreeqc.h" #include "phqalloc.h" diff --git a/phreeqcpp/class_main.cpp b/phreeqcpp/class_main.cpp index fe3d0e77..03fb3e40 100644 --- a/phreeqcpp/class_main.cpp +++ b/phreeqcpp/class_main.cpp @@ -323,7 +323,7 @@ write_banner(void) /* date */ #ifdef NPP - len = snprintf(buffer, sizeof(buffer), "%s", "July 5, 2021"); + len = snprintf(buffer, sizeof(buffer), "%s", "March 20, 2023"); #else len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@"); #endif @@ -507,7 +507,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, output_msg(sformatf(" Input file: %s\n", in_file.c_str())); output_msg(sformatf(" Output file: %s\n", out_file.c_str())); #ifdef NPP - output_msg(sformatf("Using PHREEQC: version 3.7.1, compiled July 5, 2021\n")); + output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 20, 2023\n")); #endif output_msg(sformatf("Database file: %s\n\n", token.c_str())); #ifdef NPP diff --git a/phreeqcpp/common/Utils.cxx b/phreeqcpp/common/Utils.cxx index ad96a147..8495cd36 100644 --- a/phreeqcpp/common/Utils.cxx +++ b/phreeqcpp/common/Utils.cxx @@ -10,7 +10,8 @@ #include "Utils.h" #include "Parser.h" #include "float.h" -#include "math.h" +//#include +#include #if defined(PHREEQCI_GUI) #ifdef _DEBUG diff --git a/phreeqcpp/gases.cpp b/phreeqcpp/gases.cpp index b373068a..72d8fe12 100644 --- a/phreeqcpp/gases.cpp +++ b/phreeqcpp/gases.cpp @@ -624,6 +624,7 @@ int Phreeqc:: calc_fixed_volume_gas_pressures(void) /* ---------------------------------------------------------------------- */ { + //int n_g = 0; LDBLE lp; class rxn_token *rxn_ptr; class phase *phase_ptr; @@ -644,6 +645,7 @@ calc_fixed_volume_gas_pressures(void) { if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0) PR = true; + //n_g++; } gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gas_unknowns[i]->moles); } diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index 9d07f324..053dbfb4 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -5,8 +5,8 @@ /* ---------------------------------------------------------------------- * #define DEFINITIONS * ---------------------------------------------------------------------- */ -#ifndef NAN -# define NAN -99999999 +#if !defined(NAN) +#define NAN nan("1") #endif #define MISSING -9999.999 #include "NA.h" /* NA = not available */ diff --git a/phreeqcpp/integrate.cpp b/phreeqcpp/integrate.cpp index 1cfab335..03e1c304 100644 --- a/phreeqcpp/integrate.cpp +++ b/phreeqcpp/integrate.cpp @@ -25,9 +25,9 @@ calc_all_g(void) if (use.Get_surface_ptr() == NULL) return (OK); -/* - * calculate g for each surface - */ + /* + * calculate g for each surface + */ epsilon = convergence_tolerance; if (convergence_tolerance >= 1e-8) { @@ -45,7 +45,7 @@ calc_all_g(void) if (x[j]->type != SURFACE_CB) continue; if (debug_diffuse_layer == TRUE) - output_msg(sformatf( "Calc_all_g, X[%d]\n", j)); + output_msg(sformatf("Calc_all_g, X[%d]\n", j)); cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); std::map temp_g_map; cxxSurfDL temp_g; @@ -56,11 +56,11 @@ calc_all_g(void) /* 1000 J/kJ and 1000 L/m**3 */ //alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 * // tk_x * 0.5); - alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 * - tk_x * 0.5); -/* - * calculate g for given surface for each species - */ + alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 * + tk_x * 0.5); + /* + * calculate g for given surface for each species + */ for (int i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type > HPLUS) @@ -70,10 +70,10 @@ calc_all_g(void) z_global = s_x[i]->z; if (charge_ptr->Get_grams() > 0.0) { - + if ((use.Get_surface_ptr()->Get_only_counter_ions() == false) || (((x[j]->master[0]->s->la > 0) && (z_global < 0)) - || ((x[j]->master[0]->s->la < 0) && (z_global > 0)))) + || ((x[j]->master[0]->s->la < 0) && (z_global > 0)))) { if (xd_global > 0.1) { @@ -180,11 +180,11 @@ calc_all_g(void) if (debug_diffuse_layer == TRUE) { output_msg(sformatf( - "\t%12f\t%12.4e\t%12.4e\t%12.4e\n", - (double) z_global, - (double) charge_ptr->Get_g_map()[z_global].Get_g(), - (double) new_g, - (double) (new_g - charge_ptr->Get_g_map()[z_global].Get_g()))); + "\t%12f\t%12.4e\t%12.4e\t%12.4e\n", + (double)z_global, + (double)charge_ptr->Get_g_map()[z_global].Get_g(), + (double)new_g, + (double)(new_g - charge_ptr->Get_g_map()[z_global].Get_g()))); } } charge_ptr->Get_g_map()[z_global].Set_g(new_g); @@ -201,7 +201,7 @@ calc_all_g(void) g_function(xd_global) / F_C_MOL; dg *= -2. / (exp(x[j]->master[0]->s->la * LOG_10) * - exp(x[j]->master[0]->s->la * LOG_10)); + exp(x[j]->master[0]->s->la * LOG_10)); if ((xd_global - 1) < 0.0) { dg *= -1.0; @@ -225,16 +225,16 @@ calc_all_g(void) if (debug_diffuse_layer == TRUE) { output_msg(sformatf("\nSurface component %d: charge,\tg,\tdg/dlny,\txd\n", - (int) charge_ptr->Get_g_map().size())); + (int)charge_ptr->Get_g_map().size())); std::map::iterator it; for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++) { output_msg(sformatf( - "\t%12f\t%12.4e\t%12.4e\t%12.4e\n", - (double) it->first, - (double) it->second.Get_g(), - (double) it->second.Get_dg(), - (double) xd_global)); + "\t%12f\t%12.4e\t%12.4e\t%12.4e\n", + (double)it->first, + (double)it->second.Get_g(), + (double)it->second.Get_dg(), + (double)xd_global)); } } } @@ -253,10 +253,10 @@ g_function(LDBLE x_value) return (0.0); sum = 0.0; ln_x_value = log(x_value); - + cxxSurfaceCharge *charge_ptr = &(use.Get_surface_ptr()->Get_surface_charges()[0]); std::map::iterator it = charge_ptr->Get_g_map().begin(); - for ( ; it != charge_ptr->Get_g_map().end(); it++) + for (; it != charge_ptr->Get_g_map().end(); it++) { it->second.Set_psi_to_z(exp(ln_x_value * it->first) - 1.0); } @@ -272,31 +272,31 @@ g_function(LDBLE x_value) sum = 0.0; sum1 = 0.0; output_msg(sformatf( - "Species\tmoles\tX**z-1\tsum\tsum charge\n")); + "Species\tmoles\tX**z-1\tsum\tsum charge\n")); for (i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type < H2O && s_x[i]->z != 0.0) { sum += s_x[i]->moles * (pow(x_value, s_x[i]->z) - 1.0); sum1 += s_x[i]->moles * s_x[i]->z; - output_msg(sformatf( "%s\t%e\t%e\t%e\t%e\n", - s_x[i]->name, (double) s_x[i]->moles, - (double) (pow((LDBLE) x_value, (LDBLE) s_x[i]->z) - - 1.0), (double) sum, (double) sum1)); + output_msg(sformatf("%s\t%e\t%e\t%e\t%e\n", + s_x[i]->name, (double)s_x[i]->moles, + (double)(pow((LDBLE)x_value, (LDBLE)s_x[i]->z) - + 1.0), (double)sum, (double)sum1)); } } - error_string = sformatf( "Negative sum in g_function, %e\t%e.", - (double) sum, (double) x_value); + error_string = sformatf("Negative sum in g_function, %e\t%e.", + (double)sum, (double)x_value); error_msg(error_string, CONTINUE); error_string = sformatf( - "Solutions must be charge balanced, charge imbalance is %e\n", - (double) sum1); + "Solutions must be charge balanced, charge imbalance is %e\n", + (double)sum1); error_msg(error_string, STOP); } return_value = (exp(ln_x_value * z_global) - - 1) / sqrt((x_value * x_value * mass_water_aq_x * sum)); + 1) / sqrt((x_value * x_value * mass_water_aq_x * sum)); return (return_value); } /* ---------------------------------------------------------------------- */ @@ -309,9 +309,9 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy) ns = 1; dif = fabs(xv - xa[1]); -/* - * Malloc work space - */ + /* + * Malloc work space + */ std::vector c, d; c.resize((size_t)n + 1); d.resize((size_t)n + 1); @@ -354,7 +354,7 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy) } *yv += *dy; -/* *yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */ + /* *yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */ } return; } @@ -376,7 +376,7 @@ midpnt(LDBLE x1, LDBLE x2, int n) { for (it = 1, j = 1; j < n - 1; j++) it *= 3; - tnm = (LDBLE) it; + tnm = (LDBLE)it; del = (x2 - x1) / (3 * tnm); ddel = del + del; xv = x1 + 0.5 * del; @@ -420,7 +420,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2) if (debug_diffuse_layer == TRUE) { output_msg(sformatf( - "Iterations in qromb_midpnt: %d\n", j)); + "Iterations in qromb_midpnt: %d\n", j)); } return (sv[j]); } @@ -436,7 +436,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2) if (debug_diffuse_layer == TRUE) { output_msg(sformatf( - "Iterations in qromb_midpnt: %d\n", j)); + "Iterations in qromb_midpnt: %d\n", j)); } return (ss); } @@ -444,7 +444,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2) } error_string = sformatf( - "\nToo many iterations integrating diffuse layer.\n"); + "\nToo many iterations integrating diffuse layer.\n"); error_msg(error_string, STOP); return (-999.9); } @@ -455,9 +455,9 @@ calc_init_g(void) { if (use.Get_surface_ptr() == NULL) return (OK); -/* - * calculate g for each surface - */ + /* + * calculate g for each surface + */ for (int j = 0; j < count_unknowns; j++) { if (x[j]->type != SURFACE_CB) @@ -468,7 +468,7 @@ calc_init_g(void) /* second 1000 is liters/m**3 */ //alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * // 1000.0 * tk_x * 0.5); - alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * + alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 * tk_x * 0.5); if (charge_ptr->Get_g_map().size() == 0) @@ -476,9 +476,9 @@ calc_init_g(void) cxxSurfDL temp_g; charge_ptr->Get_g_map()[0.0] = temp_g; } -/* - * calculate g for given surface for each species - */ + /* + * calculate g for given surface for each species + */ for (int i = 0; i < (int)this->s_x.size(); i++) { if (s_x[i]->type > HPLUS) @@ -511,7 +511,7 @@ calc_init_g(void) { int is = s_x[i]->number; - assert (is < (int) s_diff_layer.size()); + assert(is < (int)s_diff_layer.size()); // species found in diff_layer s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0); s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0); @@ -520,15 +520,15 @@ calc_init_g(void) if (debug_diffuse_layer == TRUE) { output_msg(sformatf( - "\nSurface component %d: charge,\tg,\tdg\n", - (int) charge_ptr->Get_g_map().size())); + "\nSurface component %d: charge,\tg,\tdg\n", + (int)charge_ptr->Get_g_map().size())); std::map::iterator it; for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++) { - output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n", - (double) it->first, - (double) it->second.Get_g(), - (double) it->second.Get_dg())); + output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\n", + (double)it->first, + (double)it->second.Get_g(), + (double)it->second.Get_dg())); } } } @@ -539,21 +539,21 @@ int Phreeqc:: initial_surface_water(void) /* ---------------------------------------------------------------------- */ { -/* - * In initial surface calculation, need to calculate - * mass of water in diffuse layer. - * diffuse layer water + aqueous solution water = bulk water. - * Ionic strength is fixed, so diffuse-layer water will not change - */ + /* + * In initial surface calculation, need to calculate + * mass of water in diffuse layer. + * diffuse layer water + aqueous solution water = bulk water. + * Ionic strength is fixed, so diffuse-layer water will not change + */ LDBLE debye_length, b, r, rd, ddl_limit, rd_limit, fraction, sum_surfs, l_s; LDBLE damp_aq; -/* - * Debye length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36 - * - * 1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters. - */ - //debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x) - // / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.); + /* + * Debye length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36 + * + * 1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters. + */ + //debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x) + // / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.); debye_length = (eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x) / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.); debye_length = sqrt(debye_length); @@ -561,10 +561,10 @@ initial_surface_water(void) /* ddl is at most the fraction ddl_limit of bulk water */ ddl_limit = use.Get_surface_ptr()->Get_DDL_limit(); -/* - * Loop through all surface components, calculate each H2O surface (diffuse layer), - * H2O aq, and H2O bulk (diffuse layers plus aqueous). - */ + /* + * Loop through all surface components, calculate each H2O surface (diffuse layer), + * H2O aq, and H2O bulk (diffuse layers plus aqueous). + */ if (use.Get_surface_ptr()->Get_debye_lengths() > 0) { @@ -606,14 +606,14 @@ initial_surface_water(void) mass_water_surfaces_x = use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 - ddl_limit); r = 0.002 * (mass_water_surfaces_x + - use.Get_solution_ptr()->Get_mass_water()) / sum_surfs; + use.Get_solution_ptr()->Get_mass_water()) / sum_surfs; rd_limit = (1 - sqrt(1 - ddl_limit)) * r; rd = rd_limit; use.Get_surface_ptr()->Set_thickness(rd); } else mass_water_surfaces_x = - (r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water(); + (r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water(); for (int i = 0; i < count_unknowns; i++) { if (x[i]->type != SURFACE_CB) @@ -692,14 +692,14 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr) if (use.Get_surface_ptr() == NULL) return (OK); -/* - * Find position of component in list of components - */ + /* + * Find position of component in list of components + */ -/* - * Loop through all surface components, calculate each H2O surface (diffuse layer), - * H2O aq, and H2O bulk (diffuse layers plus aqueous). - */ + /* + * Loop through all surface components, calculate each H2O surface (diffuse layer), + * H2O aq, and H2O bulk (diffuse layers plus aqueous). + */ count_elts = 0; paren_count = 0; mass_water_surface = charge_ptr->Get_mass_water(); @@ -717,9 +717,9 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr) } moles_excess = mass_water_aq_x * molality * g; moles_surface = mass_water_surface * molality + moles_excess; -/* - * Accumulate elements in diffuse layer - */ + /* + * Accumulate elements in diffuse layer + */ add_elt_list(s_x[j]->next_elt, moles_surface); } add_elt_list(s_h2o->next_elt, mass_water_surface / gfw_water); @@ -731,22 +731,33 @@ int Phreeqc:: calc_all_donnan(void) /* ---------------------------------------------------------------------- */ { - bool converge; + bool converge; int cd_m; - LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot; + LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot, co_ion; LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1; if (use.Get_surface_ptr() == NULL) return (OK); - //f_sinh = sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * - // tk_x * mu_x); f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * - tk_x * mu_x); -/* - * calculate g for each surface... - */ - if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths()) // DL_pitz + tk_x * mu_x); + bool only_count = use.Get_surface_ptr()->Get_only_counter_ions(); + bool correct_GC = use.Get_surface_ptr()->Get_correct_GC(); + /* calculate g for each surface... + */ + if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths() || + correct_GC) // DL_pitz && correct_GC initial_surface_water(); + LDBLE nDbl = 1; + if (correct_GC) + { + if ((nDbl = use.Get_surface_ptr()->Get_debye_lengths()) == 0) + { + LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3); + nDbl = use.Get_surface_ptr()->Get_thickness() / debye_length; + //use.Get_surface_ptr()->Set_debye_lengths(nDbl); + } + } + converge = TRUE; for (int j = 0; j < count_unknowns; j++) { @@ -755,10 +766,10 @@ calc_all_donnan(void) cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge); if (debug_diffuse_layer == TRUE) - output_msg(sformatf( "Calc_all_g, X[%d]\n", j)); -/* - * sum eq of each charge number in solution... - */ + output_msg(sformatf("Calc_all_g, X[%d]\n", j)); + /* + * sum eq of each charge number in solution... + */ std::map::iterator it; for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) { @@ -778,15 +789,18 @@ calc_all_donnan(void) f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */ f_psi = f_psi / 2; cd_m = 1; - } else + } + else { f_psi = x[j]->master[0]->s->la * LOG_10; cd_m = -1; } surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL; - if (fabs(surf_chrg_eq) > 5e3) + LDBLE lim_seq = 5e3; + if (correct_GC) lim_seq = 5e1; + if (fabs(surf_chrg_eq) > lim_seq) { - surf_chrg_eq = (surf_chrg_eq < 0 ? -5e3 : 5e3); + surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq); var1 = surf_chrg_eq / (A_surf * f_sinh / F_C_MOL); var1 = (var1 + sqrt(var1 * var1 + 1)); f_psi = (var1 > 1e-8 ? log(var1) : -18.4); @@ -799,8 +813,11 @@ calc_all_donnan(void) surf_chrg_eq2 = A_surf * f_sinh * sinh(f_psi2) / F_C_MOL; /* find psi_avg that matches surface charge... */ - psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq); - psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2); + std::vector zcorr(charge_group_map.size()); + std::vector zcorr2(charge_group_map.size()); + //LDBLE fD = 0; + psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, zcorr); + psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, zcorr2); /*output_msg(sformatf( "psi's %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */ @@ -808,9 +825,15 @@ calc_all_donnan(void) ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; ratio_aq_tot = charge_ptr->Get_mass_water() / mass_water_bulk_x; + int z_iter = 0; for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) { - LDBLE z = it->first; + LDBLE z = it->first, z1 = z; + co_ion = surf_chrg_eq * z; + if (correct_GC) + z1 = zcorr[z_iter]; + //z1 *= cgc[0] * pow(z_factor, abs(z)); + if (!ratio_aq) { charge_ptr->Get_g_map()[z].Set_g(0); @@ -819,17 +842,13 @@ calc_all_donnan(void) converge = true; continue; } - new_g = ratio_aq * (exp(cd_m * z * psi_avg) - 1); - if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0) - //((surf_chrg_eq < 0 && z < 0) - // || (surf_chrg_eq > 0 && z > 0))) + new_g = ratio_aq * (exp(cd_m * z1 * psi_avg) - 1); + if (only_count && co_ion > 0) new_g = -ratio_aq; if (new_g <= -ratio_aq) new_g = -ratio_aq + G_TOL * 1e-3; - new_g2 = ratio_aq * (exp(cd_m * z * psi_avg2) - 1); - if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0) - //((surf_chrg_eq < 0 && z < 0) - // || (surf_chrg_eq > 0 && z > 0))) + new_g2 = ratio_aq * (exp(cd_m * z1 * psi_avg2) - 1); + if (only_count && co_ion > 0) new_g2 = -ratio_aq; if (new_g2 <= -ratio_aq) new_g2 = -ratio_aq + G_TOL * 1e-3; @@ -859,9 +878,9 @@ calc_all_donnan(void) /* save Boltzmann factor * water fraction for MCD calc's in transport */ if (converge) { - if (use.Get_surface_ptr()->Get_only_counter_ions()) + if (only_count) { - if (surf_chrg_eq * z > 0) // co-ions are not in the DL + if (co_ion > 0) // co-ions are not in the DL charge_ptr->Get_z_gMCD_map()[z] = 0; else // assume that counter-ions have the free water conc for diffusion charge_ptr->Get_z_gMCD_map()[z] = ratio_aq_tot; @@ -869,23 +888,27 @@ calc_all_donnan(void) else charge_ptr->Get_z_gMCD_map()[z] = (new_g / ratio_aq + 1) * ratio_aq_tot; } + z_iter++; } if (debug_diffuse_layer == TRUE) { - std::string name = x[j]->master[0]->elt->name; + std::string name = x[j]->master[0]->elt->name; Utilities::replace("_psi", "", name); output_msg(sformatf( - "\nDonnan all on %s (%d): charge, \tg, \tdg, Psi_surface = %8f V. \n", - name.c_str(), (int) charge_ptr->Get_g_map().size(), - x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * - tk_x / F_KJ_V_EQ)); + "\nDonnan all on %s (%d): charge, \tg, \tdg, \tzcorr, \tPsi_surface = %8f V, \tDebye lengths = %8f. \n", + name.c_str(), (int)charge_ptr->Get_g_map().size(), + x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, + nDbl)); + int i = 0; for (std::map::iterator i_it = charge_ptr->Get_g_map().begin(); i_it != charge_ptr->Get_g_map().end(); i_it++) { - output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n", - (double) i_it->first, - (double) i_it->second.Get_g(), - (double) i_it->second.Get_dg())); + output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n", + (double)i_it->first, + (double)i_it->second.Get_g(), + (double)i_it->second.Get_dg(), + (double)zcorr[i])); + i++; } } } @@ -904,7 +927,7 @@ calc_init_donnan(void) //sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * // tk_x * mu_x); sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * - tk_x * mu_x); + tk_x * mu_x); if (convergence_tolerance >= 1e-8) { G_TOL = 1e-9; @@ -913,9 +936,9 @@ calc_init_donnan(void) { G_TOL = 1e-13; } -/* - * sum eq of each charge number in solution... - */ + /* + * sum eq of each charge number in solution... + */ charge_group_map.clear(); charge_group_map[0.0] = 0.0; @@ -932,9 +955,10 @@ calc_init_donnan(void) charge_group_map[s_x[i]->z] = s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl; } } -/* - * calculate g for each surface... - */ + std::vector zcorr(charge_group_map.size()); + /* + * calculate g for each surface... + */ for (int j = 0; j < count_unknowns; j++) { if (x[j]->type != SURFACE_CB) @@ -948,14 +972,13 @@ calc_init_donnan(void) { f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */ f_psi = f_psi / 2; - } else + } + else f_psi = x[j]->master[0]->s->la * LOG_10; surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL; /* find psi_avg that matches surface charge... */ -/* psi_avg = calc_psi_avg (0); - appt 7/9/8... may have to change above one */ - psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq); + psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, zcorr); /* fill in g's */ ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; @@ -975,7 +998,7 @@ calc_init_donnan(void) if (charge_ptr->Get_g_map()[z].Get_g() != 0) { - charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) / + charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) / (eq * F_C_MOL)); } else @@ -986,7 +1009,7 @@ calc_init_donnan(void) for (int i = 0; i < (int)this->s_x.size(); i++) { int is = s_x[i]->number; - assert (is < (int) s_diff_layer.size()); + assert(is < (int)s_diff_layer.size()); s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0.0); s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0.0); @@ -997,17 +1020,17 @@ calc_init_donnan(void) std::string name = x[j]->master[0]->elt->name; Utilities::replace("_psi", "", name); output_msg(sformatf( - "\nDonnan init on %s : charge, \tg, \tdg, Psi_surface = %8f V. \n", - name.c_str(), - x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * - tk_x / F_KJ_V_EQ)); + "\nDonnan init on %s : charge, \tg, \tdg, \tzcorr, Psi_surface = %8f V. \n", + name.c_str(), + x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * + tk_x / F_KJ_V_EQ)); for (std::map::iterator i_it = charge_ptr->Get_g_map().begin(); i_it != charge_ptr->Get_g_map().end(); i_it++) { - output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n", - (double) i_it->first, - (double) i_it->second.Get_g(), - (double) i_it->second.Get_dg())); + output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n", + (double)i_it->first, + (double)i_it->second.Get_g(), + (double)i_it->second.Get_dg())); } } } @@ -1015,13 +1038,13 @@ calc_init_donnan(void) } /* ---------------------------------------------------------------------- */ LDBLE Phreeqc:: -calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq) +calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector &zcorr) /* ---------------------------------------------------------------------- */ { -/* - * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge - */ - LDBLE fd, fd1, p, temp, ratio_aq; + /* + * calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge + */ + LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, z1, z1_c, eq, co_ion, sum_counter, sum_co; ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x; p = 0; @@ -1031,30 +1054,79 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq) p = -0.5 * log(-surf_chrg_eq * ratio_aq / mu_x + 1); else if (surf_chrg_eq > 0) p = 0.5 * log(surf_chrg_eq * ratio_aq / mu_x + 1); -/* - * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq - * g(p) = exp(-p * z_i) * ratio_aq - * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0: - * g(p) = (exp(-p *z_i) - 1) * ratio_aq - */ - int l_iter = 0; + /* + * Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq + * g(p) = exp(-p * z_i) * ratio_aq + * Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0: + * g(p) = (exp(-p *z_i) - 1) * ratio_aq + * with correct_GC true: + * correct ions to better match the integrated concentrations: + z == 1? z *= 0.285 cgc[6] + z == 2? z *= 0.372 cgc[7] + z == -1? z *= cgc[0] * (mu_x**( cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4]) + z == -2? z *= cgc[0] * (mu_x**(cgc[5] * cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4]) + */ + + cxxSurface *surf_p = use.Get_surface_ptr(); + bool correct_GC = surf_p->Get_correct_GC(), local_correct_GC = correct_GC; + bool only_count = surf_p->Get_only_counter_ions(); + LDBLE Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6, + cgc[10] = { 0.36, 0.1721, 0.798, 0.287, 0.1457, 1.2, 0.285, 0.372 }; + + if (!surf_p->Donnan_factors.empty()) + std::copy(std::begin(surf_p->Donnan_factors), std::end(surf_p->Donnan_factors), cgc); + + cgc[1] *= pow(nDbl, cgc[2]) * pow(Gamma, cgc[3]) * pow(mu_x, cgc[4]); + + int l_iter = 0, z_iter; + sum_co = sum_counter = 0; do { fd = surf_chrg_eq; fd1 = 0.0; + z_iter = 0; + if (l_iter == 1 && local_correct_GC && fabs(sum_counter) < fabs(sum_co)) + { + local_correct_GC = false; + l_iter = 0; + } std::map::iterator it; for (it = charge_group_map.begin(); it != charge_group_map.end(); it++) { - LDBLE z = it->first; - if (!z || (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)) + z = it->first; + z1 = z; + if (l_iter == 0) zcorr[z_iter] = z; + co_ion = surf_chrg_eq * z; + if (!z || (only_count && co_ion > 0)) + { + z_iter++; continue; - LDBLE eq = it->second; - /* multiply with ratio_aq for multiplier options cp and cm - in calc_all_donnan (not used now)... */ - temp = exp(-z * p) * ratio_aq; + } + if (nDbl && local_correct_GC) + { + /*psi_DL = fabs(p * p_psi);*/ + if (co_ion < 0) + {//counter-ion + if (fabs(z) > 1) temp = cgc[7]; + else temp = cgc[6]; + sum_counter += z * temp; + } + else + {// co-ion + if (fabs(z) > 1) temp = cgc[0] * pow(mu_x, cgc[1] * cgc[5]); + else temp = cgc[0] * pow(mu_x, cgc[1]); + sum_co += z * temp; + } + zcorr[z_iter] = z * temp; + } + z1 = zcorr[z_iter]; + eq = it->second; + temp = exp(-z1 * p) * ratio_aq; fd += eq * temp; - fd1 -= z * eq * temp; + fd1 -= z1 * eq * temp; + if (z == 1) z1_c = z1; + z_iter++; } fd /= -fd1; p += (fd > 1) ? 1 : ((fd < -1) ? -1 : fd); @@ -1076,12 +1148,11 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq) (double)surf_chrg_eq, (double)charge_ptr->Get_mass_water()); error_msg(error_string, STOP); } - } - while (fabs(fd) > 1e-12 && p != 0.0); + } while (fabs(fd) > 1e-12 && p != 0.0); if (debug_diffuse_layer == TRUE) output_msg(sformatf( - "iter in calc_psi_avg = %d. g(+1) = %8f. surface charge = %12.4e.\n", - l_iter, (double) (exp(-p) - 1), (double) surf_chrg_eq)); + "iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n", + l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * z1_c * p_psi))); return (p); } diff --git a/phreeqcpp/inverse.cpp b/phreeqcpp/inverse.cpp index b9bbcd0d..7c0b63c9 100644 --- a/phreeqcpp/inverse.cpp +++ b/phreeqcpp/inverse.cpp @@ -3586,8 +3586,10 @@ check_isotopes(class inverse *inv_ptr) if (j < inv_ptr->i_u[i].uncertainties.size() && !isnan(inv_ptr->i_u[i].uncertainties[j])) #else + //if (j < inv_ptr->i_u[i].uncertainties.size() + // && inv_ptr->i_u[i].uncertainties[j] != NAN) if (j < inv_ptr->i_u[i].uncertainties.size() - && inv_ptr->i_u[i].uncertainties[j] != NAN) + && !std::isnan(inv_ptr->i_u[i].uncertainties[j])) #endif { kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]); @@ -3598,8 +3600,10 @@ check_isotopes(class inverse *inv_ptr) else if (inv_ptr->i_u[i].uncertainties.size() > 0 && !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1])) #else + //else if (inv_ptr->i_u[i].uncertainties.size() > 0 + // && inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN) else if (inv_ptr->i_u[i].uncertainties.size() > 0 - && inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN) + && !std::isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1])) #endif { kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]); @@ -3609,7 +3613,8 @@ check_isotopes(class inverse *inv_ptr) #ifdef NPP else if (!isnan(kit->second.Get_ratio_uncertainty())) #else - else if (kit->second.Get_ratio_uncertainty() != NAN) + //else if (kit->second.Get_ratio_uncertainty() != NAN) + else if (!std::isnan(kit->second.Get_ratio_uncertainty())) #endif { kit->second.Set_x_ratio_uncertainty( @@ -3641,7 +3646,8 @@ check_isotopes(class inverse *inv_ptr) #ifdef NPP if (isnan(kit->second.Get_x_ratio_uncertainty())) #else - if (kit->second.Get_x_ratio_uncertainty() == NAN) + //if (kit->second.Get_x_ratio_uncertainty() == NAN) + if (std::isnan(kit->second.Get_x_ratio_uncertainty())) #endif { error_string = sformatf( diff --git a/phreeqcpp/isotopes.cpp b/phreeqcpp/isotopes.cpp index 68570473..46e3c318 100644 --- a/phreeqcpp/isotopes.cpp +++ b/phreeqcpp/isotopes.cpp @@ -903,7 +903,8 @@ punch_calculate_values(void) #ifdef NPP if (isnan(rate_moles)) #else - if (rate_moles == NAN) + //if (rate_moles == NAN) + if (std::isnan(rate_moles)) #endif { error_string = sformatf( "Calculated value not SAVEed for %s.", @@ -1136,7 +1137,8 @@ calculate_values(void) #ifdef NPP if (isnan(rate_moles)) #else - if (rate_moles == NAN) + //if (rate_moles == NAN) + if (std::isnan(rate_moles)) #endif { error_string = sformatf( "Calculated value not SAVEed for %s.", @@ -1203,7 +1205,8 @@ calculate_values(void) #ifdef NPP if (isnan(rate_moles)) #else - if (rate_moles == NAN) + //if (rate_moles == NAN) + if (std::isnan(rate_moles)) #endif { error_string = sformatf( "Calculated value not SAVEed for %s.", diff --git a/phreeqcpp/kinetics.cpp b/phreeqcpp/kinetics.cpp index f6e02338..6e59bcb1 100644 --- a/phreeqcpp/kinetics.cpp +++ b/phreeqcpp/kinetics.cpp @@ -118,7 +118,8 @@ calc_kinetic_reaction(cxxKinetics *kinetics_ptr, LDBLE time_step) #ifdef NPP if (isnan(rate_moles)) #else - if (rate_moles == NAN) + //if (rate_moles == NAN) + if (std::isnan(rate_moles)) #endif { error_string = sformatf( "Moles of reaction not SAVEed for %s.", diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 5cc76a3f..baa20e71 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -540,6 +540,8 @@ initial_exchangers(int print) viscosity(); species_list_sort(); print_exchange(); + if (pr.user_print) + print_user_print(); xexchange_save(n_user); punch_all(); /* free_model_allocs(); */ @@ -671,7 +673,9 @@ initial_gas_phases(int print) } print_gas_phase(); - if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/) + if (pr.user_print) + print_user_print(); + if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/) warning_msg("While initializing gas phase composition by equilibrating:\n" " Found definitions of gas` critical temperature and pressure.\n" " Going to use Peng-Robinson in subsequent calculations.\n"); @@ -745,7 +749,8 @@ initial_surfaces(int print) set_and_run_wrapper(-1, FALSE, FALSE, -1, 0.0); species_list_sort(); print_surface(); - /*print_all(); */ + if (pr.user_print) + print_user_print(); punch_all(); xsurface_save(n_user); /* free_model_allocs(); */ @@ -1257,6 +1262,7 @@ xsolution_save(int n_user) temp_solution.Set_ah2o(ah2o_x); //temp_solution.Set_density(density_x); temp_solution.Set_density(calc_dens()); + temp_solution.Set_viscosity(this->viscos); temp_solution.Set_total_h(total_h_x); temp_solution.Set_total_o(total_o_x); temp_solution.Set_cb(cb_x); /* cb_x does not include surface charge sfter sum_species */ diff --git a/phreeqcpp/model.cpp b/phreeqcpp/model.cpp index bf8d8fd7..eb5be6ab 100644 --- a/phreeqcpp/model.cpp +++ b/phreeqcpp/model.cpp @@ -2543,6 +2543,7 @@ int Phreeqc:: calc_gas_pressures(void) /* ---------------------------------------------------------------------- */ { + //int n_g = 0; LDBLE lp, V_m = 0; class rxn_token *rxn_ptr; std::vector phase_ptrs; @@ -2575,6 +2576,7 @@ calc_gas_pressures(void) phase_ptrs.push_back(phase_ptr); if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0) PR = true; + //n_g++; } if (iterations > 2 && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME) { @@ -3776,6 +3778,7 @@ residuals(void) int converge; LDBLE l_toler; + //LDBLE sum_residual; LDBLE sinh_constant; LDBLE sum, sum1; class master *master_ptr, *master_ptr1, *master_ptr2; @@ -3783,6 +3786,7 @@ residuals(void) int print_fail; std::vector cd_psi; print_fail = FALSE; + //sum_residual = 0.0; sigmaddl = 0; sum = 0; /* @@ -4526,6 +4530,7 @@ residuals(void) * Store residuals in array */ my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; + //sum_residual += fabs(residual[i]); } /* * Return diff --git a/phreeqcpp/pitzer.cpp b/phreeqcpp/pitzer.cpp index b71f24d3..35d494fc 100644 --- a/phreeqcpp/pitzer.cpp +++ b/phreeqcpp/pitzer.cpp @@ -52,6 +52,7 @@ pitzer_tidy(void) const char *string1, *string2; int i, j, order; int i0, i1, i2; + //int count_pos, count_neg, count_neut, count[3], jj; int count_neut, count[3], jj; LDBLE z0, z1; class pitz_param *pzp_ptr; @@ -339,21 +340,22 @@ pitzer_tidy(void) i0 = pitz_params[i]->ispec[0]; i1 = pitz_params[i]->ispec[1]; i2 = pitz_params[i]->ispec[2]; + //count_pos = count_neg = count_neut = 0; count_neut = 0; for (j = 0; j <= 2; j++) { - //if (spec[pitz_params[i]->ispec[j]]->z > 0) - //{ - // count_pos++; - //} + if (spec[pitz_params[i]->ispec[j]]->z > 0) + { + //count_pos++; + } if (spec[pitz_params[i]->ispec[j]]->z == 0) { count_neut++; } - //if (spec[pitz_params[i]->ispec[j]]->z < 0) - //{ - // count_neg++; - //} + if (spec[pitz_params[i]->ispec[j]]->z < 0) + { + //count_neg++; + } } /* All neutral */ if (count_neut == 3) diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp index 0d26ba98..f98d584c 100644 --- a/phreeqcpp/prep.cpp +++ b/phreeqcpp/prep.cpp @@ -1980,6 +1980,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr) * Output: space is allocated and a list of master species pointers is * returned. */ + //int j, l, count_list; int j, l; char token[MAX_LENGTH]; std::vector master_ptr_list; @@ -1987,6 +1988,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr) /* * Make list of master species pointers */ + //count_list = 0; //master_ptr_list = unknown_alloc_master(); master_ptr0 = master_ptr; if (master_ptr0 == master_ptr->s->primary) @@ -2146,6 +2148,7 @@ mb_for_species_aq(int n) * by coef, usually moles. * mb_unknowns.coef - coefficient of s[n] in equation or relation */ + //int i, j; int i; class master *master_ptr; class unknown *unknown_ptr; @@ -2222,6 +2225,7 @@ mb_for_species_aq(int n) */ if (use.Get_surface_ptr() != NULL && s[n]->type < H2O && dl_type_x != cxxSurface::NO_DL) { + //j = 0; for (i = 0; i < count_unknowns; i++) { if (x[i]->type == SURFACE_CB) @@ -2233,6 +2237,7 @@ mb_for_species_aq(int n) store_mb_unknowns(unknown_ptr, s_diff_layer[n][charge_ptr->Get_name()].Get_g_moles_address(), s[n]->z, s_diff_layer[n][charge_ptr->Get_name()].Get_dg_g_moles_address()); + //j++; } } } diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp index 612460be..e32714e9 100644 --- a/phreeqcpp/print.cpp +++ b/phreeqcpp/print.cpp @@ -339,8 +339,12 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr) { LDBLE exp_g = charge_ptr->Get_g_map()[1].Get_g() * mass_water_aq_x / mass_water_surface + 1; LDBLE psi_DL = -log(exp_g) * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ; - output_msg(sformatf( - "\n\tTotal moles in diffuse layer (excluding water), Donnan calculation.")); + if (use.Get_surface_ptr()->Get_correct_GC()) + output_msg(sformatf( + "\n\tTotal moles in diffuse layer (excluding water), Donnan corrected to match Poisson-Boltzmann.")); + else + output_msg(sformatf( + "\n\tTotal moles in diffuse layer (excluding water), Donnan calculation.")); output_msg(sformatf( "\n\tDonnan Layer potential, psi_DL = %10.3e V.\n\tBoltzmann factor, exp(-psi_DL * F / RT) = %9.3e (= c_DL / c_free if z is +1).\n\n", psi_DL, exp_g)); @@ -2234,7 +2238,7 @@ print_totals(void) (double) calc_solution_volume())); } /* VP: Density End */ -#ifdef NPP +//#ifdef NPP if (print_viscosity) { output_msg(sformatf("%45s%9.5f", "Viscosity (mPa s) = ", @@ -2251,7 +2255,7 @@ print_totals(void) } else output_msg(sformatf("\n")); } -#endif +//#endif output_msg(sformatf("%45s%7.3f\n", "Activity of water = ", exp(s_h2o->la * LOG_10))); output_msg(sformatf("%45s%11.3e\n", "Ionic strength (mol/kgw) = ", diff --git a/phreeqcpp/read.cpp b/phreeqcpp/read.cpp index ba2a0b1c..3052877d 100644 --- a/phreeqcpp/read.cpp +++ b/phreeqcpp/read.cpp @@ -6283,7 +6283,7 @@ read_surface(void) * ERROR if error occurred reading data * */ - int n_user; + int n_user, i1; LDBLE conc; const char* cptr, *cptr1; std::string token, token1, name; @@ -6306,9 +6306,10 @@ read_surface(void) "ccm", /* 13 */ "equilibrium", /* 14 */ "site_units", /* 15 */ - "ddl" /* 16 */ + "ddl", /* 16 */ + "donnan_factors" /* 17 */ }; - int count_opt_list = 17; + int count_opt_list = 18; /* * kin_surf is for Surfaces, related to kinetically reacting minerals * they are defined if "sites" is followed by mineral name: @@ -6413,7 +6414,7 @@ read_surface(void) if (thickness != 0) { error_msg - ("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8", + ("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5 limit 0.9 correct_GC true\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8 correct_GC false", CONTINUE); error_msg(line_save, CONTINUE); input_error++; @@ -6436,6 +6437,23 @@ read_surface(void) break; } } + else if (token[0] == 'C' || token[0] == 'c') + { + copy_token(token1, &next_char); + if (token1[0] == 'T' || token1[0] == 't' || token1[0] == 'F' || token1[0] == 'f') + { + temp_surface.Set_correct_GC(get_true_false(token1.c_str(), TRUE) == TRUE); + continue; + } else + { + error_msg + ("Expected True or False for correct_GC (which brings co-ion concentrations closer to their integrated double layer value).", + CONTINUE); + error_msg(line_save, CONTINUE); + input_error++; + break; + } + } else if (token[0] == 'V' || token[0] == 'v') { int j = copy_token(token1, &next_char); @@ -6570,6 +6588,31 @@ read_surface(void) case 16: /* ddl */ temp_surface.Set_type(cxxSurface::DDL); break; + case 17: /* Donnan_factors */ + temp_surface.Donnan_factors.clear(); + i1 = 0; + for (;;) + { + int i = copy_token(token, &next_char); + if (i == DIGIT && i1 < 8) + { + (void)sscanf(token.c_str(), SCANFORMAT, &dummy); + temp_surface.Donnan_factors.push_back(dummy); + i1++; + continue; + } + else if (i != EMPTY || i1 > 8) + { + error_msg + ("Expected at most 8 numbers for the Donnan_factors for co- and counter-ions,\n z *= cgc[0] * (mu_x**(cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * mu_x**(cgc[4])", + CONTINUE); + error_msg(line_save, CONTINUE); + input_error++; + break; + } + break; + } + break; case OPTION_DEFAULT: /* * Read surface component diff --git a/phreeqcpp/spread.cpp b/phreeqcpp/spread.cpp index 5d7137b5..1f442138 100644 --- a/phreeqcpp/spread.cpp +++ b/phreeqcpp/spread.cpp @@ -1141,6 +1141,7 @@ copy_token_tab(std::string& token, const char **cptr) * EOL, * UNKNOWN. */ + //int i, return_value; int return_value; char c; /* @@ -1180,6 +1181,7 @@ copy_token_tab(std::string& token, const char **cptr) /* * Begin copying to token */ + //i = 0; for (;;) { c = **cptr; @@ -1196,6 +1198,7 @@ copy_token_tab(std::string& token, const char **cptr) { token.push_back(c); (*cptr)++; + //i++; } } return (return_value); diff --git a/phreeqcpp/sundialsmath.cpp b/phreeqcpp/sundialsmath.cpp index 0f0c4578..42487bb7 100644 --- a/phreeqcpp/sundialsmath.cpp +++ b/phreeqcpp/sundialsmath.cpp @@ -59,7 +59,8 @@ #include -#include + //#include +#include #include "sundialsmath.h" #include "sundialstypes.h" diff --git a/phreeqcpp/tally.cpp b/phreeqcpp/tally.cpp index 131411e1..8fdb0e6a 100644 --- a/phreeqcpp/tally.cpp +++ b/phreeqcpp/tally.cpp @@ -798,6 +798,7 @@ build_tally_table(void) */ int j, k, l, p, save_print_use; size_t n; + //int count_tt_pure_phase, count_tt_ss_phase, count_tt_kinetics; class phase *phase_ptr; char token[MAX_LENGTH]; const char* cptr; @@ -870,6 +871,7 @@ build_tally_table(void) /* * Count pure phases */ + //count_tt_pure_phase = 0; if (Rxn_pp_assemblage_map.size() > 0) { /* @@ -902,6 +904,7 @@ build_tally_table(void) /* * Add to table */ + //count_tt_pure_phase++; n = count_tally_table_columns; extend_tally_table(); tally_table[n].name = phase_ptr->name; @@ -928,6 +931,7 @@ build_tally_table(void) /* * Add solid-solution pure phases */ + //count_tt_ss_phase = 0; if (Rxn_ss_assemblage_map.size() > 0) { /* @@ -960,6 +964,7 @@ build_tally_table(void) /* * Add to table */ + //count_tt_ss_phase++; n = count_tally_table_columns; extend_tally_table(); tally_table[n].name = phase_ptr->name; @@ -977,6 +982,7 @@ build_tally_table(void) /* * Add kinetic reactants */ + //count_tt_kinetics = 0; if (Rxn_kinetics_map.size() > 0) { std::map::iterator it; @@ -1000,6 +1006,7 @@ build_tally_table(void) /* * Add to table */ + //count_tt_kinetics++; n = count_tally_table_columns; extend_tally_table(); tally_table[n].name = string_hsave(kinetics_comp_ptr->Get_rate_name().c_str()); diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp index 71280107..db5574b1 100644 --- a/phreeqcpp/tidy.cpp +++ b/phreeqcpp/tidy.cpp @@ -982,7 +982,8 @@ tidy_gas_phase(void) #ifdef NPP if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) #else - if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN) + //if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN) + if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) #endif { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); @@ -1015,7 +1016,8 @@ tidy_gas_phase(void) #ifdef NPP if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) #else - if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN) + //if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN) + if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read())) #endif { P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read(); @@ -1690,7 +1692,8 @@ tidy_ss_assemblage(void) #ifdef NPP if (isnan(comp_ptr->Get_moles())) #else - if (comp_ptr->Get_moles() == NAN) + //if (comp_ptr->Get_moles() == NAN) + if (std::isnan(comp_ptr->Get_moles())) #endif { input_error++; @@ -3021,7 +3024,8 @@ tidy_isotopes(void) #ifdef NPP if (!isnan(master[k]->isotope_ratio_uncertainty)) #else - if (master[k]->isotope_ratio_uncertainty != NAN) + //if (master[k]->isotope_ratio_uncertainty != NAN) + if (!std::isnan(master[k]->isotope_ratio_uncertainty)) #endif { temp_iso.Set_ratio_uncertainty_defined(true); diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 36dda3cc..378d818b 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -1888,16 +1888,17 @@ fill_spec(int l_cell_no, int ref_cell) if (por_il < interlayer_Dpor_lim) por_il = viscos_il_f = 0.0; /* - * correct diffusion coefficient for temperature and viscosity, D_T = D_298 * Tk * viscos_298 / (298 * viscos) + * correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos * modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959 */ viscos = viscos_0; /* * put temperature factor in por_factor which corrects for porous medium... */ - viscos_f *= tk_x * viscos_0_25 / (298.15 * viscos); - viscos_il_f *= tk_x * viscos_0_25 / (298.15 * viscos); - sol_D[l_cell_no].viscos_f = tk_x * viscos_0_25 / (298.15 * viscos); + dum = viscos_0_25 / viscos; + viscos_f *= dum; + viscos_il_f *= dum; + sol_D[l_cell_no].viscos_f = dum; count_spec = count_exch_spec = 0; /* @@ -1908,6 +1909,8 @@ fill_spec(int l_cell_no, int ref_cell) qsort(&species_list[0], species_list.size(), sizeof(class species_list), sort_species_name); } + if (correct_Dw) + calc_SC(); for (i = 0; i < (int)species_list.size(); i++) { @@ -2097,7 +2100,7 @@ fill_spec(int l_cell_no, int ref_cell) } if (correct_Dw) { - calc_SC(); // note that neutral species are corrected as if z = 1, but is viscosity-dependent + //calc_SC(); // removed that neutral species are corrected as if z = 1, but is viscosity-dependent sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr * viscos_f; } if (l_cell_no <= count_cells + 1 && sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn) > diffc_max) @@ -2962,20 +2965,20 @@ diffuse_implicit(LDBLE DDt, int stagnant) if (!strcmp(ct[icell].m_s[cp].name, "H")) { dummy = ct[icell].m_s[cp].tot1; - if (dV_dcell || (icell > 0 && icell <= ilast)) + if (dV_dcell || fix_current || (icell > 0 && icell <= ilast)) sptr1->Set_total_h(sptr1->Get_total_h() - dummy); - if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) sptr2->Set_total_h(sptr2->Get_total_h() + dummy); if (sptr_stag) { dummy = ct[icell].m_s[cp].tot_stag; - if (dV_dcell || (icell > 0 && icell <= ilast)) + if (dV_dcell || fix_current || (icell > 0 && icell <= ilast)) sptr1->Set_total_h(sptr1->Get_total_h() - dummy); if (icell == c) { // mix the boundary solution with the stagnant column end-cell dummy += ct[icell + 1].m_s[cp].tot_stag; - if (dV_dcell) + if (dV_dcell || fix_current) sptr2->Set_total_h(sptr2->Get_total_h() - ct[icell + 1].m_s[cp].tot_stag); } sptr_stag->Set_total_h(sptr_stag->Get_total_h() + dummy); @@ -2985,19 +2988,19 @@ diffuse_implicit(LDBLE DDt, int stagnant) if (!strcmp(ct[icell].m_s[cp].name, "O")) { dummy = ct[icell].m_s[cp].tot1; - if (dV_dcell || (icell > 0 && icell <= ilast)) + if (dV_dcell || fix_current || (icell > 0 && icell <= ilast)) sptr1->Set_total_o(sptr1->Get_total_o() - dummy); - if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) sptr2->Set_total_o(sptr2->Get_total_o() + dummy); if (sptr_stag) { dummy = ct[icell].m_s[cp].tot_stag; - if (dV_dcell || (icell > 0 && icell <= ilast)) + if (dV_dcell || fix_current || (icell > 0 && icell <= ilast)) sptr1->Set_total_o(sptr1->Get_total_o() - dummy); if (icell == c) { dummy += ct[icell + 1].m_s[cp].tot_stag; - if (dV_dcell) + if (dV_dcell || fix_current) sptr2->Set_total_o(sptr2->Get_total_o() - ct[icell + 1].m_s[cp].tot_stag); } sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy); @@ -3018,7 +3021,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) // check for negative moles, add moles from other redox states and the donnan layer when necessary and available... if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol && - (dV_dcell || (icell > 0 && icell <= ilast))) + (dV_dcell || fix_current || (icell > 0 && icell <= ilast))) { dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name); if (ct[icell].dl_s > 1e-8) @@ -3039,7 +3042,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) if (icell == c && sptr_stag && ct[c1].m_s[cp].tot_stag) dum = ct[c1].m_s[cp].tot_stag; if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol && - (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))) + (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))) { dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name); if (ct[icell + 1].dl_s > 1e-8) @@ -3057,7 +3060,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) if (fabs(ct[icell].m_s[cp].tot2) < fabs(ct[icell].m_s[cp].tot1)) ct[icell].m_s[cp].tot1 = ct[icell].m_s[cp].tot2; - if (dV_dcell || (icell > 0 && icell <= ilast)) + if (dV_dcell || fix_current || (icell > 0 && icell <= ilast)) { dum = ct[icell].m_s[cp].tot1; if (stagnant) @@ -3078,10 +3081,10 @@ diffuse_implicit(LDBLE DDt, int stagnant) //ct[icell].J_ij_sum -= dum * ct[icell].m_s[cp].charge; } - if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) { dum = ct[icell].m_s[cp].tot1; - if (stagnant && icell == c && dV_dcell) + if (stagnant && icell == c && (dV_dcell || fix_current)) dum -= ct[c1].m_s[cp].tot_stag; dum2 += dum; sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol); @@ -3134,7 +3137,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) } // reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted... - if (dV_dcell && icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1) + if ((dV_dcell/* || fix_current*/) && icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1) { dummy = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_totals()[ct[0].m_s[cp].name] / ct[0].kgw * (1 - ct[0].dl_s); if (dummy > 1e-6) @@ -4828,10 +4831,10 @@ Step (from cell 1 to count_cells + 1): cxxSurface *surface_ptr1, *surface_ptr2; LDBLE viscos_f; /* - * temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos) + * temperature and viscosity correction for MCD coefficient, D_T = D_298 * viscos_298 / viscos */ viscos_f = viscos_0; - viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f); + viscos_f = viscos_0_25 / viscos_f; //n1 = 0; //n2 = n1 + 1; @@ -5502,7 +5505,7 @@ diff_stag_surf(int mobile_cell) * temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos) */ viscos_f = viscos_0; - viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f); + viscos_f = viscos_0_25 / viscos_f; cxxSurface surface_n1, surface_n2; cxxSurface *surface_n1_ptr = &surface_n1; @@ -5991,8 +5994,7 @@ viscosity(void) } // find B * m and D * m * mu^d3 Bc += (s_x[i]->Jones_Dole[0] + - s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * - t1; + s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1; // define f_I from the exponent of the D * m^d3 term... if (s_x[i]->Jones_Dole[5] >= 1) t2 = mu_x / 3 / s_x[i]->Jones_Dole[5]; @@ -6074,8 +6076,8 @@ viscosity(void) viscos += (viscos_0 * (Bc + Dc)); if (viscos < 0) { - viscos = 0; - warning_msg("viscosity < 0, reset to 0."); + viscos = viscos_0; + warning_msg("viscosity < 0, reset to viscosity of pure water"); } return viscos;