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 '7939dd32097aa6395b7bd6b8e335b710a49f2b4a'

Browse Source
This commit is contained in:
github-actions[bot] 2025-09-10 22:46:08 +00:00
commit eca261b232
6 changed files with 126 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/phreeqcpp/NameDouble.cxx
View File

@ -226,7 +226,7 @@ cxxNameDouble::add_intensive(const cxxNameDouble & addee, LDBLE f1,
// Sums two name doubles, this*f1 + f2*nd2
//
{
assert(f1 >= 0 && f2 >= 0);
//assert(f1 >= 0 && f2 >= 0);
for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end();
it++)
{
@ -248,7 +248,7 @@ cxxNameDouble::add_log_activities(const cxxNameDouble & addee, LDBLE f1,
// Sums two name doubles, this*f1 + f2*nd2, assuming log values
//
{
assert(f1 >= 0 && f2 >= 0);
//assert(f1 >= 0 && f2 >= 0);
for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end();
it++)
{

43
src/phreeqcpp/Solution.cxx
View File

@ -130,23 +130,32 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
this->potV = cxxsoln_ptr1->potV;
}
//
// Sort to enable positive mixes first
const std::map < int, LDBLE >& mixcomps = mix.Get_mixComps();
std::set<std::pair<LDBLE, int> >s;
for (std::map < int, LDBLE >::const_iterator it = mixcomps.begin();
it != mixcomps.end(); it++)
{
std::pair<LDBLE, int> p(it->second, it->first);
s.insert(p);
}
//
// Mix solutions
//
const std::map < int, LDBLE >&mixcomps = mix.Get_mixComps();
std::map < int, LDBLE >::const_iterator it;
for (it = mixcomps.begin(); it != mixcomps.end(); it++)
std::set < std::pair< LDBLE, int > >::const_reverse_iterator rit = s.rbegin();
for (rit = s.rbegin(); rit != s.rend(); rit++)
{
sol = solutions.find(it->first);
sol = solutions.find(rit->second);
if (sol == solutions.end())
{
std::ostringstream msg;
msg << "Solution " << it->first << " not found in mix_cxxSolutions.";
msg << "Solution " << rit->second << " not found in mix_cxxSolutions.";
error_msg(msg.str(), CONTINUE);
}
else
{
cxxsoln_ptr1 = &(sol->second);
this->add(*cxxsoln_ptr1, it->second);
this->add(*cxxsoln_ptr1, rit->first);
}
}
@ -1418,8 +1427,25 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
return;
LDBLE ext1 = this->mass_water;
LDBLE ext2 = addee.mass_water * extensive;
LDBLE f1 = ext1 / (ext1 + ext2);
LDBLE f2 = ext2 / (ext1 + ext2);
this->mass_water += addee.mass_water * extensive;
if (this->mass_water <= 0.0)
{
std::ostringstream msg;
msg << "Negative mass of water when mixing solutions.";
error_msg(msg.str(), STOP);
}
LDBLE fconc = 0.0, f1 = 0.0, f2 = 0.0;
if (extensive > 0.0)
{
f1 = ext1 / (ext1 + ext2);
f2 = ext2 / (ext1 + ext2);
}
else
{
f1 = 1.0;
f2 = 0.0;
}
this->tc = f1 * this->tc + f2 * addee.tc;
this->ph = f1 * this->ph + f2 * addee.ph;
this->pe = f1 * this->pe + f2 * addee.pe;
@ -1433,7 +1459,6 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
this->viscos_0 = f1 * this->viscos_0 + f2 * addee.viscos_0;
this->patm = f1 * this->patm + f2 * addee.patm;
// this->potV = f1 * this->potV + f2 * addee.potV; // appt
this->mass_water += addee.mass_water * extensive;
this->soln_vol += addee.soln_vol * extensive;
this->total_alkalinity += addee.total_alkalinity * extensive;
this->totals.add_extensive(addee.totals, extensive);

3
src/phreeqcpp/global_structures.h
View File

@ -1586,7 +1586,7 @@ public:
typedef enum
{ TYPE_B0, TYPE_B1, TYPE_B2, TYPE_C0, TYPE_THETA, TYPE_LAMBDA, TYPE_ZETA,
TYPE_PSI, TYPE_ETHETA, TYPE_ALPHAS, TYPE_MU, TYPE_ETA, TYPE_Other,
TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI
TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI, TYPE_SIT_EPSILON2
} pitz_param_type;
class pitz_param
{
@ -1625,6 +1625,7 @@ public:
LDBLE eps;
LDBLE eps1;
LDBLE aphi;
LDBLE eps2;
} U;
LDBLE a[6];
LDBLE alpha;

2
src/phreeqcpp/pitzer.cpp
View File

@ -2690,7 +2690,7 @@ pitzer_make_lists(void)
/*
TYPE_B0, TYPE_B1, TYPE_B2, TYPE_C0, TYPE_THETA, TYPE_LAMBDA, TYPE_ZETA,
TYPE_PSI, TYPE_ETHETA, TYPE_ALPHAS, TYPE_MU, TYPE_ETA, TYPE_Other,
TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU
TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI, TYPE_SIT_EPSILON2
*/
int i0 = pitz_params[i]->ispec[0];
int i1 = pitz_params[i]->ispec[1];

29
src/phreeqcpp/sit.cpp
View File

@ -188,9 +188,10 @@ read_sit(void)
const char* next_char;
const char *opt_list[] = {
"epsilon", /* 0 */
"epsilon1" /* 1 */
"epsilon1", /* 1 */
"epsilon2" /* 2 */
};
int count_opt_list = 2;
int count_opt_list = 3;
/*
* Read lines
*/
@ -237,6 +238,11 @@ read_sit(void)
n = 2;
opt_save = OPTION_DEFAULT;
break;
case 2: /* epsilon2 */
pzp_type = TYPE_SIT_EPSILON2;
n = 2;
opt_save = OPTION_DEFAULT;
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
@ -275,6 +281,9 @@ calc_sit_param(class pitz_param *pz_ptr, LDBLE TK, LDBLE TR)
case TYPE_SIT_EPSILON_MU:
pz_ptr->U.eps1 = param;
break;
case TYPE_SIT_EPSILON2:
pz_ptr->U.eps2 = param;
break;
case TYPE_Other:
default:
error_msg("Should not be TYPE_Other in function calc_sit_param",
@ -399,6 +408,7 @@ sit(void)
* Sums for sit_LGAMMA, and OSMOT
* epsilons are tabulated for log10 gamma (not ln gamma)
*/
LDBLE logmu = log10(I);
for (size_t j = 0; j < param_list.size(); j++)
{
int i = param_list[j];
@ -425,7 +435,7 @@ sit(void)
case TYPE_SIT_EPSILON_MU:
sit_LGAMMA[i0] += sit_M[i1] * I * param;
sit_LGAMMA[i1] += sit_M[i0] * I * param;
OSMOT += sit_M[i0] * sit_M[i1] * param;
//OSMOT += sit_M[i0] * sit_M[i1] * param;
if (z0 == 0.0 && z1 == 0.0)
{
OSMOT += sit_M[i0] * sit_M[i1] * param * I / 2.0;
@ -435,6 +445,19 @@ sit(void)
OSMOT += sit_M[i0] * sit_M[i1] * param * I;
}
break;
case TYPE_SIT_EPSILON2:
sit_LGAMMA[i0] += sit_M[i1] * logmu * param;
sit_LGAMMA[i1] += sit_M[i0] * logmu * param;
//OSMOT += sit_M[i0] * sit_M[i1] * param;
if (z0 == 0.0 && z1 == 0.0)
{
OSMOT += sit_M[i0] * sit_M[i1] * param * logmu / 2.0;
}
else
{
OSMOT += sit_M[i0] * sit_M[i1] * param * logmu;
}
break;
default:
case TYPE_Other:
error_msg("TYPE_Other in pitz_param list.", STOP);

120
src/phreeqcpp/step.cpp
View File

@ -621,70 +621,72 @@ add_mix(cxxMix *mix_ptr)
* calls add_solution to accumulate all data in master->totals
* and other variables.
*/
LDBLE sum_fractions, intensive, extensive;
LDBLE sum_fractions_water=0, sum_positive_water=0, intensive_water=0, extensive_water=0;
cxxSolution *solution_ptr;
int count_positive;
LDBLE sum_positive;
//LDBLE sum_fractions, intensive, extensive;
//LDBLE sum_fractions_water=0, sum_positive_water=0, intensive_water=0, extensive_water=0;
//cxxSolution *solution_ptr;
//int count_positive;
//LDBLE sum_positive;
if (mix_ptr == NULL)
return (OK);
if (mix_ptr->Get_mixComps().size() == 0)
return (OK);
sum_fractions = 0.0;
sum_positive = 0.0;
count_positive = 0;
std::map<int, LDBLE>::const_iterator it;
for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
{
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
if (solution_ptr == NULL)
{
error_string = sformatf("Mix solution not found, %d.",
it->first);
error_msg(error_string, CONTINUE);
input_error++;
continue;
}
sum_fractions += it->second;
sum_fractions_water += it->second * solution_ptr->Get_mass_water();
if (it->second > 0)
{
sum_positive += it->second;
sum_positive_water += it->second * solution_ptr->Get_mass_water();
count_positive++;
}
}
for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
{
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
if (solution_ptr == NULL)
{
error_string = sformatf( "Mix solution not found, %d.",
it->first);
error_msg(error_string, CONTINUE);
input_error++;
continue;
}
extensive = it->second;
extensive_water = it->second * solution_ptr->Get_mass_water();
intensive = extensive / sum_fractions;
intensive_water = extensive_water / sum_fractions_water;
if (count_positive < (int) mix_ptr->Get_mixComps().size())
{
if (it->second > 0)
{
intensive = extensive / sum_positive;
intensive_water = extensive_water / sum_positive_water;
}
else
{
intensive = 0;
}
}
add_solution(solution_ptr, extensive, intensive_water);
}
//sum_fractions = 0.0;
//sum_positive = 0.0;
//count_positive = 0;
//std::map<int, LDBLE>::const_iterator it;
//for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
//{
// solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
// if (solution_ptr == NULL)
// {
// error_string = sformatf("Mix solution not found, %d.",
// it->first);
// error_msg(error_string, CONTINUE);
// input_error++;
// continue;
// }
// sum_fractions += it->second;
// sum_fractions_water += it->second * solution_ptr->Get_mass_water();
// if (it->second > 0)
// {
// sum_positive += it->second;
// sum_positive_water += it->second * solution_ptr->Get_mass_water();
// count_positive++;
// }
//}
//for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
//{
// solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
// if (solution_ptr == NULL)
// {
// error_string = sformatf( "Mix solution not found, %d.",
// it->first);
// error_msg(error_string, CONTINUE);
// input_error++;
// continue;
// }
// extensive = it->second;
// extensive_water = it->second * solution_ptr->Get_mass_water();
// intensive = extensive / sum_fractions;
// intensive_water = extensive_water / sum_fractions_water;
//
// if (count_positive < (int) mix_ptr->Get_mixComps().size())
// {
// if (it->second > 0)
// {
// intensive = extensive / sum_positive;
// intensive_water = extensive_water / sum_positive_water;
// }
// else
// {
// intensive = 0;
// }
// }
// add_solution(solution_ptr, extensive, intensive_water);
//}
cxxSolution s(Rxn_solution_map, *mix_ptr, 1);
add_solution(&s, 1.0, 1.0);
return (OK);
}
#ifdef SKIP_ERROR