max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-15 16:18:22 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge commit 'b34eedb91d08cd3b5fd34a70311fad8c14c7854d'

Browse Source
This commit is contained in:
Darth Vader 2024-11-14 01:37:56 +00:00
commit 57b3a6cd39
8 changed files with 229 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
phreeqcpp/Phreeqc.cpp
View File

@ -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;

3
phreeqcpp/Phreeqc.h
View File

@ -473,6 +473,7 @@ public:
#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 @@ public:
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<int, cxxMix>& mix_map);
//int read_solution_mix(void);
@ -1674,6 +1676,7 @@ protected:
std::vector<double> x_arg, res_arg, scratch;
/* gases.cpp ------------------------------- */
LDBLE a_aa_sum, b2, b_sum, R_TK;
std::map < std::pair<std::string, std::string>, double > gas_binary_parameters;
/* input.cpp ------------------------------- */
int check_line_return;

4
phreeqcpp/PhreeqcKeywords/Keywords.cpp
View File

@ -134,6 +134,7 @@ std::map<const std::string, Keywords::KEYWORDS>::value_type("rate_parameters_pk"
std::map<const std::string, Keywords::KEYWORDS>::value_type("rate_parameters_svd", Keywords::KEY_RATE_PARAMETERS_SVD),
std::map<const std::string, Keywords::KEYWORDS>::value_type("rate_parameters_hermanska", Keywords::KEY_RATE_PARAMETERS_HERMANSKA),
std::map<const std::string, Keywords::KEYWORDS>::value_type("mean_gammas", Keywords::KEY_MEAN_GAMMAS),
std::map<const std::string, Keywords::KEYWORDS>::value_type("gas_binary_parameters", Keywords::KEY_GAS_BINARY_PARAMETERS),
std::map<const std::string, Keywords::KEYWORDS>::value_type("solution_mix", Keywords::KEY_SOLUTION_MIX),
std::map<const std::string, Keywords::KEYWORDS>::value_type("mix_solution", Keywords::KEY_SOLUTION_MIX),
std::map<const std::string, Keywords::KEYWORDS>::value_type("exchange_mix", Keywords::KEY_EXCHANGE_MIX),
@ -228,7 +229,8 @@ std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_REACTI
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_RATE_PARAMETERS_PK, "RATE_PARAMETERS_PK"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_RATE_PARAMETERS_SVD, "RATE_PARAMETERS_SVD"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_RATE_PARAMETERS_HERMANSKA, "RATE_PARAMETERS_HERMANSKA"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_MEAN_GAMMAS, "RATE_MEAN_GAMMAS"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_MEAN_GAMMAS, "MEAN_GAMMAS"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_GAS_BINARY_PARAMETERS, "GAS_BINARY_PARAMETERS"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_SOLUTION_MIX, "SOLUTION_MIX"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_EXCHANGE_MIX, "EXCHANGE_MIX"),
std::map<Keywords::KEYWORDS, const std::string>::value_type(Keywords::KEY_GAS_PHASE_MIX, "GAS_PHASE_MIX"),

1
phreeqcpp/PhreeqcKeywords/Keywords.h
View File

@ -80,6 +80,7 @@ public:
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,

6
phreeqcpp/basicsubs.cpp
View File

@ -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);

116
phreeqcpp/gases.cpp
View File

@ -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<std::pair<std::string, std::string>, 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;
}

64
phreeqcpp/prep.cpp
View File

@ -3843,36 +3843,37 @@ calc_PR(std::vector<class phase *> 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<class phase *> 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;

100
phreeqcpp/read.cpp
View File

@ -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<std::string, std::string> 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)
/* ---------------------------------------------------------------------- */
{