From b2885f29f98cc04fc19ae3e130fa9187ea696fc1 Mon Sep 17 00:00:00 2001 From: Darth Vader Date: Thu, 14 Nov 2024 01:42:25 +0000 Subject: [PATCH] Squashed 'src/' changes from ccb9ba3c..f5587dae f5587dae Made GAS_BINARY_PARAMETERS backward compatible f57db333 GAS_BINARY_PARAMETERS implemented. gas_binary_parameters test case f4f1f698 Tweaked H+ viscosity in databases, new test cases fig1 and H_HCl_HBr git-subtree-dir: src git-subtree-split: f5587dae0a068b1346a04108be71e266860e4706 --- Phreeqc.cpp | 1 + Phreeqc.h | 3 + PhreeqcKeywords/Keywords.cpp | 4 +- PhreeqcKeywords/Keywords.h | 1 + basicsubs.cpp | 6 +- gases.cpp | 116 +++++++++++++++++++++++++---------- prep.cpp | 64 +++++++++---------- read.cpp | 100 ++++++++++++++++++++++++++++++ 8 files changed, 229 insertions(+), 66 deletions(-) diff --git a/Phreeqc.cpp b/Phreeqc.cpp index 672b22d2c..c5c892c4b 100644 --- a/Phreeqc.cpp +++ b/Phreeqc.cpp @@ -1756,6 +1756,7 @@ Phreeqc::InternalCopy(const Phreeqc* pSrc) b2 = pSrc->b2; b_sum = pSrc->b_sum; R_TK = pSrc->R_TK; + gas_binary_parameters = pSrc->gas_binary_parameters; /* input.cpp ------------------------------- */ check_line_return = 0; reading_db = FALSE; diff --git a/Phreeqc.h b/Phreeqc.h index ddff8262b..cd5425168 100644 --- a/Phreeqc.h +++ b/Phreeqc.h @@ -473,6 +473,7 @@ class Phreeqc #endif int calc_gas_pressures(void); int calc_fixed_volume_gas_pressures(void); + double calc_gas_binary_parameter(std::string name1, std::string name2) const; int calc_ss_fractions(void); int gammas(LDBLE mu); int gammas_a_f(int i); @@ -700,6 +701,7 @@ class Phreeqc int read_rate_parameters_svd(void); int read_rate_parameters_hermanska(void); int read_mean_gammas(void); + int read_gas_binary_parameters(void); int read_mix(void); int read_entity_mix(std::map& mix_map); //int read_solution_mix(void); @@ -1674,6 +1676,7 @@ class Phreeqc std::vector x_arg, res_arg, scratch; /* gases.cpp ------------------------------- */ LDBLE a_aa_sum, b2, b_sum, R_TK; + std::map < std::pair, double > gas_binary_parameters; /* input.cpp ------------------------------- */ int check_line_return; diff --git a/PhreeqcKeywords/Keywords.cpp b/PhreeqcKeywords/Keywords.cpp index 8b0e8c000..590fc597c 100644 --- a/PhreeqcKeywords/Keywords.cpp +++ b/PhreeqcKeywords/Keywords.cpp @@ -134,6 +134,7 @@ std::map::value_type("rate_parameters_pk" std::map::value_type("rate_parameters_svd", Keywords::KEY_RATE_PARAMETERS_SVD), std::map::value_type("rate_parameters_hermanska", Keywords::KEY_RATE_PARAMETERS_HERMANSKA), std::map::value_type("mean_gammas", Keywords::KEY_MEAN_GAMMAS), +std::map::value_type("gas_binary_parameters", Keywords::KEY_GAS_BINARY_PARAMETERS), std::map::value_type("solution_mix", Keywords::KEY_SOLUTION_MIX), std::map::value_type("mix_solution", Keywords::KEY_SOLUTION_MIX), std::map::value_type("exchange_mix", Keywords::KEY_EXCHANGE_MIX), @@ -228,7 +229,8 @@ std::map::value_type(Keywords::KEY_REACTI std::map::value_type(Keywords::KEY_RATE_PARAMETERS_PK, "RATE_PARAMETERS_PK"), std::map::value_type(Keywords::KEY_RATE_PARAMETERS_SVD, "RATE_PARAMETERS_SVD"), std::map::value_type(Keywords::KEY_RATE_PARAMETERS_HERMANSKA, "RATE_PARAMETERS_HERMANSKA"), -std::map::value_type(Keywords::KEY_MEAN_GAMMAS, "RATE_MEAN_GAMMAS"), +std::map::value_type(Keywords::KEY_MEAN_GAMMAS, "MEAN_GAMMAS"), +std::map::value_type(Keywords::KEY_GAS_BINARY_PARAMETERS, "GAS_BINARY_PARAMETERS"), std::map::value_type(Keywords::KEY_SOLUTION_MIX, "SOLUTION_MIX"), std::map::value_type(Keywords::KEY_EXCHANGE_MIX, "EXCHANGE_MIX"), std::map::value_type(Keywords::KEY_GAS_PHASE_MIX, "GAS_PHASE_MIX"), diff --git a/PhreeqcKeywords/Keywords.h b/PhreeqcKeywords/Keywords.h index 06de5596e..ae10e6fa0 100644 --- a/PhreeqcKeywords/Keywords.h +++ b/PhreeqcKeywords/Keywords.h @@ -80,6 +80,7 @@ class Keywords KEY_RATE_PARAMETERS_SVD, KEY_RATE_PARAMETERS_HERMANSKA, KEY_MEAN_GAMMAS, + KEY_GAS_BINARY_PARAMETERS, KEY_SOLUTION_MIX, KEY_EXCHANGE_MIX, KEY_GAS_PHASE_MIX, diff --git a/basicsubs.cpp b/basicsubs.cpp index 917353dac..adac27361 100644 --- a/basicsubs.cpp +++ b/basicsubs.cpp @@ -420,7 +420,7 @@ calc_SC(void) q = 1 / ((t1 / z_plus + (1 - t1) / z_min) * (z_min + z_plus)); sqrt_q = sqrt(q); - // B1 = relaxtion, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka)) + // B1 = relaxation, B2 = electrophoresis in ll = (ll0 - B2 * sqrt(mu) / f2(1 + ka)) * (1 - B1 * sqrt(mu) / f1(1 + ka)) a = 1.60218e-19 * 1.60218e-19 / (6 * pi); B1 = a / (2 * 8.8542e-12 * eps_r * 1.38066e-23 * tk_x) * q / (1 + sqrt_q) * DH_B * 1e10 * z_plus * z_min; // DH_B is per Angstrom (*1e10) B2 = a * AVOGADRO / viscos_0 * DH_B * 1e17; // DH_B per Angstrom (*1e10), viscos in mPa.s (*1e3), B2 in cm2 (*1e4) @@ -486,8 +486,10 @@ calc_SC(void) //av += 0 * t1; } Dw *= Dw_SC * l_z; - if (!a2 || !strcmp(s_x[i]->name, "H+")) + if (!a2) t1 = 1; + else if (!strcmp(s_x[i]->name, "H+")) + t1 = pow(1 + mu_x, a2); else { v0 = calc_vm0(s_x[i]->name, tc_x, 1, 0); diff --git a/gases.cpp b/gases.cpp index bfcc64d74..dea1e551d 100644 --- a/gases.cpp +++ b/gases.cpp @@ -432,36 +432,37 @@ calc_PR(void) continue; a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * phase_ptr1->pr_a * phase_ptr1->pr_alpha); - if (!strcmp(phase_ptr->name, "H2O(g)")) - { - if (!strcmp(phase_ptr1->name, "CO2(g)")) - a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Propane(g)")) - a_aa *= 0.45; - } - if (!strcmp(phase_ptr1->name, "H2O(g)")) - { - if (!strcmp(phase_ptr->name, "CO2(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Propane(g)")) - a_aa *= 0.45; - } + a_aa *= calc_gas_binary_parameter(phase_ptr->name, phase_ptr1->name); + //if (!strcmp(phase_ptr->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr1->name, "CO2(g)")) + // a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + // else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Propane(g)")) + // a_aa *= 0.45; + //} + //if (!strcmp(phase_ptr1->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr->name, "CO2(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Propane(g)")) + // a_aa *= 0.45; + //} a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; a_aa_sum2 += phase_ptr1->fraction_x * a_aa; } @@ -566,12 +567,14 @@ calc_PR(void) if (ri + rq / 2 <= 0) { V_m = pow(ri - rq / 2, one_3) + pow(- ri - rq / 2, one_3) - r3[1] / 3; - } else + } + else { ri = - pow(ri + rq / 2, one_3); V_m = ri - rp / (3.0 * ri) - r3[1] / 3; } - } else // use complex plane... + } + else // use complex plane... { ri = sqrt(- rp3 / 27); // rp < 0 ri1 = acos(- rq / 2 / ri); @@ -695,3 +698,52 @@ calc_fixed_volume_gas_pressures(void) return (OK); } +/* ---------------------------------------------------------------------- */ +double Phreeqc:: +calc_gas_binary_parameter(std::string name1, std::string name2) const +/* ---------------------------------------------------------------------- */ +{ + double f = 1.0; + std::pair < std::string, std::string > p; + p = { name1, name2 }; + std::map, double>::const_iterator gas_pair_it; + gas_pair_it = gas_binary_parameters.find(p); + if (gas_pair_it != gas_binary_parameters.end()) + { + f = (1.0 - gas_pair_it->second); + } + else + { + if (!strcmp(name1.c_str(), "H2O(g)")) + { + if (!strcmp(name2.c_str(), "CO2(g)")) + f = 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + else if (!strcmp(name2.c_str(), "H2S(g)") || !strcmp(name2.c_str(), "H2Sg(g)")) + f = 0.81; + else if (!strcmp(name2.c_str(), "CH4(g)") || !strcmp(name2.c_str(), "Mtg(g)") || !strcmp(name2.c_str(), "Methane(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "N2(g)") || !strcmp(name2.c_str(), "Ntg(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "Ethane(g)")) + f = 0.51; + else if (!strcmp(name2.c_str(), "Propane(g)")) + f = 0.45; + } + if (!strcmp(name2.c_str(), "H2O(g)")) + { + if (!strcmp(name1.c_str(), "CO2(g)")) + f = 0.81; + else if (!strcmp(name1.c_str(), "H2S(g)") || !strcmp(name1.c_str(), "H2Sg(g)")) + f = 0.81; + else if (!strcmp(name1.c_str(), "CH4(g)") || !strcmp(name1.c_str(), "Mtg(g)") || !strcmp(name1.c_str(), "Methane(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "N2(g)") || !strcmp(name1.c_str(), "Ntg(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "Ethane(g)")) + f = 0.51; + else if (!strcmp(name1.c_str(), "Propane(g)")) + f = 0.45; + } + } + return f; +} \ No newline at end of file diff --git a/prep.cpp b/prep.cpp index fb867e087..9ab2668e9 100644 --- a/prep.cpp +++ b/prep.cpp @@ -3843,36 +3843,37 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) continue; a_aa = sqrt(phase_ptr->pr_a * phase_ptr->pr_alpha * phase_ptr1->pr_a * phase_ptr1->pr_alpha); - if (!strcmp(phase_ptr->name, "H2O(g)")) - { - if (!strcmp(phase_ptr1->name, "CO2(g)")) - a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Propane(g)")) - a_aa *= 0.45; - } - if (!strcmp(phase_ptr1->name, "H2O(g)")) - { - if (!strcmp(phase_ptr->name, "CO2(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) - a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Ethane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Propane(g)")) - a_aa *= 0.45; - } + a_aa *= calc_gas_binary_parameter(phase_ptr->name, phase_ptr1->name); + //if (!strcmp(phase_ptr->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr1->name, "CO2(g)")) + // a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 + // else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr1->name, "Propane(g)")) + // a_aa *= 0.45; + //} + //if (!strcmp(phase_ptr1->name, "H2O(g)")) + //{ + // if (!strcmp(phase_ptr->name, "CO2(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) + // a_aa *= 0.81; + // else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Ethane(g)")) + // a_aa *= 0.51; + // else if (!strcmp(phase_ptr->name, "Propane(g)")) + // a_aa *= 0.45; + //} a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; a_aa_sum2 += phase_ptr1->fraction_x * a_aa; } @@ -3946,7 +3947,8 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) } if (P <= 0) // iterations = -1 P = 1; - } else + } + else { if (P < 1e-10) P = 1e-10; diff --git a/read.cpp b/read.cpp index bb451db68..1b0bb84d8 100644 --- a/read.cpp +++ b/read.cpp @@ -147,6 +147,9 @@ read_input(void) case Keywords::KEY_MEAN_GAMMAS: read_mean_gammas(); break; + case Keywords::KEY_GAS_BINARY_PARAMETERS: + read_gas_binary_parameters(); + break; case Keywords::KEY_SOLUTION_MIX: //read_solution_mix(); read_entity_mix(Rxn_solution_mix_map); @@ -2551,6 +2554,103 @@ read_mean_gammas(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +read_gas_binary_parameters(void) +/* ---------------------------------------------------------------------- */ +{ + /* + * Reads GAS_BINARY_PARAMETERS data + * + * Arguments: + * none + * + * Returns: + * KEYWORD if keyword encountered, input_error may be incremented if + * a keyword is encountered in an unexpected position + * EOF if eof encountered while reading mass balance concentrations + * ERROR if error occurred reading data + * + */ + std::string token; + int return_value, opt; + const char* next_char; + const char* opt_list[] = { + "xxxx", /* 0 */ + }; + int count_opt_list = 0; + /* + * Read rate parameters + */ + return_value = UNKNOWN; + for (;;) + { + opt = get_option(opt_list, count_opt_list, &next_char); + switch (opt) + { + case OPTION_EOF: /* end of file */ + return_value = EOF; + break; + case OPTION_KEYWORD: /* keyword */ + return_value = KEYWORD; + break; + case OPTION_DEFAULT: /* add to gas_binary_parameters map */ + { + bool error = false; + std::string gas1, gas2; + int j = copy_token(token, &next_char); + if (j != EMPTY) + { + gas1 = token; + } + else + { + error = true; + } + j = copy_token(token, &next_char); + if (j != EMPTY) + { + gas2 = token; + } + else + { + error = true; + } + j = copy_token(token, &next_char); + double d; + if (j != EMPTY) + { + j = sscanf(token.c_str(), SCANFORMAT, &d); + } + else + { + error = true; + } + if (!error) + { + std::pair p; + p = { gas1, gas2 }; + gas_binary_parameters[p] = d; + p = { gas2, gas1 }; + gas_binary_parameters[p] = d; + } + else + { + error_msg("Error reading gas binary parameter", CONTINUE); + } + } + break; + case OPTION_ERROR: + input_error++; + error_msg("Unknown input in GAS_BINARY_PARAMETERS keyword.", CONTINUE); + error_msg(line_save, CONTINUE); + break; + } + if (return_value == EOF || return_value == KEYWORD) + break; + } + return (return_value); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: read_rate_parameters_svd(void) /* ---------------------------------------------------------------------- */ {