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 remote-tracking branch 'origin/master' into gtest

Browse Source
This commit is contained in:
Charlton, Scott R 2020-09-07 22:25:59 -06:00
commit b030a0b9db
10 changed files with 899 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
phreeqcpp/ExchComp.cxx
View File

@ -406,7 +406,7 @@ cxxExchComp::multiply(LDBLE extensive)
{
this->totals.multiply(extensive);
this->charge_balance *= extensive;
this->phase_proportion *= extensive;
//this->phase_proportion *= extensive;
}
void

13
phreeqcpp/PBasic.cpp
View File

@ -1612,6 +1612,9 @@ listtokens(FILE * f, tokenrec * l_buf)
case tokdh_a:
output_msg("DH_A"); // Debye-Hueckel A
break;
case tokdebye_length:
output_msg("DEBYE_LENGTH"); // Debye-Hueckel length
break;
case tokdh_b:
output_msg("DH_B"); // Debye-Hueckel B
break;
@ -3706,6 +3709,13 @@ factor(struct LOC_exec * LINK)
n.UU.val = PhreeqcPtr->DH_A;
}
break;
case tokdebye_length:
{
double debye_length = (PhreeqcPtr->eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * PhreeqcPtr->tk_x)
/ (2. * F_C_MOL * F_C_MOL * PhreeqcPtr->mu_x * 1000.);
n.UU.val = sqrt(debye_length);
break;
}
case tokdh_b:
if (PhreeqcPtr->llnl_count_temp > 0)
{
@ -3993,7 +4003,7 @@ factor(struct LOC_exec * LINK)
// call callback Basic function
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->basic_callback(x1, x2, str);
PhreeqcPtr->PHRQ_free(str);
}
break;
@ -7377,6 +7387,7 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("eps_r", PBasic::tokeps_r),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("vm", PBasic::tokvm),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_a", PBasic::tokdh_a),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("debye_length", PBasic::tokdebye_length),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_b", PBasic::tokdh_b),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_av", PBasic::tokdh_av),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("qbrn", PBasic::tokqbrn),

1
phreeqcpp/PBasic.h
View File

@ -317,6 +317,7 @@ public:
tokvm,
tokphase_vm,
tokdh_a,
tokdebye_length,
tokdh_b,
tokdh_av,
tokqbrn,

5
phreeqcpp/Phreeqc.h
View File

@ -1033,15 +1033,19 @@ public:
int tidy_logk(void);
int tidy_exchange(void);
int tidy_min_exchange(void);
int update_min_exchange(void);
int tidy_kin_exchange(void);
int update_kin_exchange(void);
int tidy_gas_phase(void);
int tidy_inverse(void);
int tidy_isotopes(void);
int tidy_isotope_ratios(void);
int tidy_isotope_alphas(void);
int tidy_kin_surface(void);
int update_kin_surface(void);
int tidy_master_isotope(void);
int tidy_min_surface(void);
int update_min_surface(void);
int tidy_phases(void);
int tidy_pp_assemblage(void);
int tidy_solutions(void);
@ -1066,6 +1070,7 @@ public:
LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);
void calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant);
void diffuse_implicit(LDBLE DDt, int stagnant);
int fill_spec(int cell_no, int ref_cell);
LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name);

30
phreeqcpp/class_main.cpp
View File

@ -109,7 +109,10 @@ main_method(int argc, char *argv[])
{
return errors;
}
#ifndef NO_UTF8_ENCODING
#ifdef DOS
SetConsoleOutputCP(CP_UTF8);
#endif
write_banner();
#endif
@ -199,7 +202,10 @@ main_method(int argc, char *argv[])
{
return errors;
}
#ifndef NO_UTF8_ENCODING
#ifdef DOS
SetConsoleOutputCP(CP_UTF8);
#endif
write_banner();
#endif
@ -271,9 +277,9 @@ write_banner(void)
char buffer[80];
int len, indent;
screen_msg(
" ÛßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÛ\n");
" █▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█\n");
screen_msg(
" º º\n");
" ║ ║\n");
/* version */
#ifdef NPP
@ -282,21 +288,21 @@ write_banner(void)
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
indent = (44 - len) / 2;
screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer,
screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer,
44 - indent - len, ' '));
screen_msg(
" º º\n");
" ║ ║\n");
screen_msg(
" º A hydrogeochemical transport model º\n");
" ║ A hydrogeochemical transport model ║\n");
screen_msg(
" º º\n");
" ║ ║\n");
screen_msg(
" º by º\n");
" ║ by ║\n");
screen_msg(
" º D.L. Parkhurst and C.A.J. Appelo º\n");
" ║ D.L. Parkhurst and C.A.J. Appelo ║\n");
screen_msg(
" º º\n");
" ║ ║\n");
/* date */
@ -306,11 +312,11 @@ write_banner(void)
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
indent = (44 - len) / 2;
screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer,
screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer,
44 - indent - len, ' '));
screen_msg(
" ÛÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÛ\n\n");
" █▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄█\n\n");
return 0;
}
@ -485,7 +491,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
}
local_database_file->close();
delete local_database_file;
user_database = (char *) free_check_null(user_database);
user_database = string_duplicate(token);
screen_msg(sformatf("Database file: %s\n\n", token));

9
phreeqcpp/mainsubs.cpp
View File

@ -2679,6 +2679,11 @@ save_init(int i)
void
Phreeqc::do_mixes(void)
{
bool surf, exch, kin, min;
surf = (Rxn_surface_mix_map.size() > 0);
exch = (Rxn_exchange_mix_map.size() > 0);
kin = (Rxn_kinetics_mix_map.size() > 0);
min = (Rxn_pp_assemblage_mix_map.size() > 0);
Utilities::Rxn_mix(Rxn_solution_mix_map, Rxn_solution_map, this);
Utilities::Rxn_mix(Rxn_exchange_mix_map, Rxn_exchange_map, this);
Utilities::Rxn_mix(Rxn_gas_phase_mix_map, Rxn_gas_phase_map, this);
@ -2686,4 +2691,8 @@ Phreeqc::do_mixes(void)
Utilities::Rxn_mix(Rxn_pp_assemblage_mix_map, Rxn_pp_assemblage_map, this);
Utilities::Rxn_mix(Rxn_ss_assemblage_mix_map, Rxn_ss_assemblage_map, this);
Utilities::Rxn_mix(Rxn_surface_mix_map, Rxn_surface_map, this);
if (exch || kin) update_kin_exchange();
if (exch || min) update_min_exchange();
if (surf || min) update_min_surface();
if (surf || kin) update_kin_surface();
}

45
phreeqcpp/print.cpp
View File

@ -1479,11 +1479,12 @@ print_species(void)
{
output_msg(sformatf("%50s%10s%10s%10s\n", "Log", "Log", "Log", "mole V"));
}
output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%10s\n\n", "Species",
#ifdef NO_UTF8_ENCODING
output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%10s\n\n", "Species",
"Molality", "Activity", "Molality", "Activity", "Gamma", "cm3/mol"));
#else
"Molality", "Activity", "Molality", "Activity", "Gamma", "cm³/mol"));
output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%11s\n\n", "Species",
"Molality", "Activity", "Molality", "Activity", "Gamma", "cm³/mol"));
#endif
/*
* Print list of species
@ -1651,7 +1652,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("\t%11.3e sigma, C/m2\n",
#else
output_msg(sformatf("\t%11.3e sigma, C/m²\n",
output_msg(sformatf("\t%11.3e sigma, C/m²\n",
#endif
(double) (charge * F_C_MOL /
(charge_ptr->Get_specific_area() *
@ -1662,7 +1663,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("\tundefined sigma, C/m2\n"));
#else
output_msg(sformatf("\tundefined sigma, C/m²\n"));
output_msg(sformatf("\tundefined sigma, C/m²\n"));
#endif
}
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
@ -1684,7 +1685,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e specific area, m2/mol %s\n",
#else
"\t%11.3e specific area, m²/mol %s\n",
"\t%11.3e specific area, m²/mol %s\n",
#endif
(double) charge_ptr->Get_specific_area(),
comp_ptr->Get_phase_name().c_str()));
@ -1692,7 +1693,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e m2 for %11.3e moles of %s\n\n",
#else
"\t%11.3e m² for %11.3e moles of %s\n\n",
"\t%11.3e m² for %11.3e moles of %s\n\n",
#endif
(double) (charge_ptr->Get_grams() *
charge_ptr->Get_specific_area()),
@ -1705,7 +1706,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e specific area, m2/mol %s\n",
#else
"\t%11.3e specific area, m²/mol %s\n",
"\t%11.3e specific area, m²/mol %s\n",
#endif
(double) charge_ptr->Get_specific_area(),
comp_ptr->Get_rate_name().c_str()));
@ -1713,7 +1714,7 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e m2 for %11.3e moles of %s\n\n",
#else
"\t%11.3e m² for %11.3e moles of %s\n\n",
"\t%11.3e m² for %11.3e moles of %s\n\n",
#endif
(double) (charge_ptr->Get_grams() *
charge_ptr->Get_specific_area()),
@ -1726,13 +1727,13 @@ print_surface(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e specific area, m2/g\n",
#else
"\t%11.3e specific area, m²/g\n",
"\t%11.3e specific area, m²/g\n",
#endif
(double) charge_ptr->Get_specific_area()));
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("\t%11.3e m2 for %11.3e g\n\n",
#else
output_msg(sformatf("\t%11.3e m² for %11.3e g\n\n",
output_msg(sformatf("\t%11.3e m² for %11.3e g\n\n",
#endif
(double) (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams()),
@ -1948,28 +1949,28 @@ print_surface_cd_music(void)
#ifdef NO_UTF8_ENCODING
"\t%11.3e sigma, plane 0, C/m2\n",
#else
"\t%11.3e sigma, plane 0, C/m²\n",
"\t%11.3e sigma, plane 0, C/m²\n",
#endif
(double) charge_ptr->Get_sigma0()));
output_msg(sformatf(
#ifdef NO_UTF8_ENCODING
"\t%11.3e sigma, plane 1, C/m2\n",
#else
"\t%11.3e sigma, plane 1, C/m²\n",
"\t%11.3e sigma, plane 1, C/m²\n",
#endif
(double) charge_ptr->Get_sigma1()));
output_msg(sformatf(
#ifdef NO_UTF8_ENCODING
"\t%11.3e sigma, plane 2, C/m2\n",
#else
"\t%11.3e sigma, plane 2, C/m²\n",
"\t%11.3e sigma, plane 2, C/m²\n",
#endif
(double) charge_ptr->Get_sigma2()));
output_msg(sformatf(
#ifdef NO_UTF8_ENCODING
"\t%11.3e sigma, diffuse layer, C/m2\n\n",
#else
"\t%11.3e sigma, diffuse layer, C/m²\n\n",
"\t%11.3e sigma, diffuse layer, C/m²\n\n",
#endif
(double) charge_ptr->Get_sigmaddl()));
}
@ -1978,7 +1979,7 @@ print_surface_cd_music(void)
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("\tundefined sigma, C/m2\n"));
#else
output_msg(sformatf("\tundefined sigma, C/m²\n"));
output_msg(sformatf("\tundefined sigma, C/m²\n"));
#endif
}
output_msg(sformatf("\t%11.3e psi, plane 0, V\n",
@ -2232,11 +2233,12 @@ print_totals(void)
if (SC > 0)
{
//output_msg(sformatf("%36s%i%7s%i\n",
output_msg(sformatf("%35s%3.0f%7s%i\n",
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("%35s%3.0f%7s%i\n",
"Specific Conductance (uS/cm, ", tc_x, "oC) = ", (int) SC));
#else
"Specific Conductance (µS/cm, ", tc_x, "°C) = ", (int) SC));
output_msg(sformatf("%36s%3.0f%7s%i\n",
"Specific Conductance (µS/cm, ", tc_x, "°C) = ", (int) SC));
#endif
}
/* VP: Density Start */
@ -2246,7 +2248,7 @@ print_totals(void)
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("%45s%9.5f", "Density (g/cm3) = ",
#else
output_msg(sformatf("%45s%9.5f", "Density (g/cm³) = ",
output_msg(sformatf("%46s%9.5f", "Density (g/cm³) = ",
#endif
(double) dens));
if (state == INITIAL_SOLUTION && use.Get_solution_ptr()->Get_initial_data()->Get_calc_density())
@ -2266,11 +2268,12 @@ print_totals(void)
(double) viscos));
if (tc_x > 200 && !pure_water)
{
output_msg(sformatf("%18s\n",
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("%18s\n",
" (solute contributions limited to 200 oC)"));
#else
" (solute contributions limited to 200 °C)"));
output_msg(sformatf("%19s\n",
" (solute contributions limited to 200 °C)"));
#endif
}
else output_msg(sformatf("\n"));
@ -2300,7 +2303,7 @@ print_totals(void)
#ifdef NO_UTF8_ENCODING
output_msg(sformatf("%45s%6.2f\n", "Temperature (oC) = ",
#else
output_msg(sformatf("%45s%6.2f\n", "Temperature (°C) = ",
output_msg(sformatf("%46s%6.2f\n", "Temperature (°C) = ",
#endif
(double) tc_x));

6
phreeqcpp/readtr.cpp
View File

@ -915,9 +915,9 @@ read_transport(void)
warning_msg(error_string);
for (i = count_por; i < all_cells - st; i++)
{
if (i == max_cells)
continue;
assert((i+1) < all_cells);
//if (i == max_cells)
// continue;
//assert((i+1) < all_cells);
if ((i+1) < all_cells)
{
cell_data[i + 1].por = pors[count_por - 1];

655
phreeqcpp/tidy.cpp
View File

@ -292,7 +292,45 @@ tidy_model(void)
{
tidy_solutions();
}
/*
* need to update exchange and surface related in case anything has changed
*/
if (keycount[Keywords::KEY_KINETICS] > 0 ||
keycount[Keywords::KEY_KINETICS_RAW] > 0 ||
keycount[Keywords::KEY_KINETICS_MODIFY] ||
keycount[Keywords::KEY_EXCHANGE] > 0 ||
keycount[Keywords::KEY_EXCHANGE_RAW] > 0 ||
keycount[Keywords::KEY_EXCHANGE_MODIFY])
{
update_kin_exchange();
}
if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] ||
keycount[Keywords::KEY_EXCHANGE] > 0 ||
keycount[Keywords::KEY_EXCHANGE_RAW] > 0 ||
keycount[Keywords::KEY_EXCHANGE_MODIFY])
{
update_min_exchange();
}
if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] ||
keycount[Keywords::KEY_SURFACE] > 0 ||
keycount[Keywords::KEY_SURFACE_RAW] > 0 ||
keycount[Keywords::KEY_SURFACE_MODIFY] > 0)
{
update_min_surface();
}
if (keycount[Keywords::KEY_KINETICS] > 0 ||
keycount[Keywords::KEY_KINETICS_RAW] > 0 ||
keycount[Keywords::KEY_KINETICS_MODIFY] > 0 ||
keycount[Keywords::KEY_SURFACE] > 0 ||
keycount[Keywords::KEY_SURFACE_RAW] > 0 ||
keycount[Keywords::KEY_SURFACE_MODIFY] > 0)
{
update_kin_surface();
}
/* if (new_model || new_exchange || new_pp_assemblage || new_surface || new_gas_phase || new_kinetics) reset_last_model(); */
if (new_model)
{
@ -3575,10 +3613,10 @@ tidy_kin_exchange(void)
assert(false);
}
cxxExchange * exchange_ptr = &(it->second);
//if (!exchange_ptr->Get_new_def())
// continue;
//if (exchange_ptr->Get_n_user() < 0)
// continue;
if (!exchange_ptr->Get_new_def())
continue;
if (exchange_ptr->Get_n_user() < 0)
continue;
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
{
@ -3669,6 +3707,140 @@ tidy_kin_exchange(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_kin_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion. Exchange needs to be updated if the
* amount of kinetic reaction has changed. Corner case of
* zero moles.
*/
{
cxxKinetics* kinetics_ptr;
char* ptr;
LDBLE conc;
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
for ( ; it != Rxn_exchange_map.end(); it++)
{
cxxExchange* exchange_ptr = &(it->second);
if (exchange_ptr->Get_n_user() < 0) continue;
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
{
cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j];
if (comp_ref.Get_rate_name().size() == 0) continue;
double comp_moles = 0.0;
/* First find exchange master species */
cxxNameDouble nd = comp_ref.Get_totals();
cxxNameDouble::iterator kit = nd.begin();
bool found_exchange = false;
for (; kit != nd.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr == NULL || elt_ptr->master == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
kit->first.c_str());
error_msg(error_string, CONTINUE);
continue;
}
if (elt_ptr->master->type == EX)
{
comp_moles = kit->second;
found_exchange = true;
}
}
//if (!found_exchange)
//{
// input_error++;
// error_string = sformatf(
// "Exchange formula does not contain an exchange master species, %s",
// comp_ref.Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
/* Now find associated kinetic reaction ... */
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, exchange_ptr->Get_n_user())) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use exchange related to kinetic reaction, %s",
exchange_ptr->Get_n_user(), comp_ref.Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
size_t k;
for (k = 0; k < kinetics_ptr->Get_kinetics_comps().size(); k++)
{
if (strcmp_nocase
(comp_ref.Get_rate_name().c_str(),
kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to exchanger, %s, not found in KINETICS %d",
comp_ref.Get_rate_name().c_str(), comp_ref.Get_formula().c_str(), exchange_ptr->Get_n_user());
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
comp_ref.Set_rate_name(kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str());
/* make exchanger concentration proportional to mineral ... */
conc = kinetics_ptr->Get_kinetics_comps()[k].Get_m() * comp_ref.Get_phase_proportion();
if (found_exchange && comp_moles > 0.0)
{
/* parse formula */
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, 1.0);
free_check_null(temp_formula);
}
cxxNameDouble nd_formula = elt_list_NameDouble();
double comp_coef = 0;
for (kit = nd_formula.begin(); kit != nd_formula.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr->master->type == EX)
{
comp_coef = kit->second;
}
}
comp_ref.multiply(comp_coef * conc / comp_moles);
}
else /* need to generate totals from scratch */
{
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comp_ref.Set_totals(elt_list_NameDouble());
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_min_exchange(void)
/* ---------------------------------------------------------------------- */
/*
@ -3694,10 +3866,10 @@ tidy_min_exchange(void)
assert(false);
}
cxxExchange * exchange_ptr = &(it->second);
//if (!exchange_ptr->Get_new_def())
// continue;
//if (exchange_ptr->Get_n_user() < 0)
// continue;
if (!exchange_ptr->Get_new_def())
continue;
if (exchange_ptr->Get_n_user() < 0)
continue;
n = exchange_ptr->Get_n_user();
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
@ -3832,6 +4004,187 @@ tidy_min_exchange(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_min_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion. Need to check in case exchange or min
* are modified.
*/
{
int n, jj;
char* ptr;
LDBLE conc;
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
for ( ; it != Rxn_exchange_map.end(); it++)
{
cxxExchange* exchange_ptr = &(it->second);
if (exchange_ptr->Get_n_user() < 0) continue;
n = exchange_ptr->Get_n_user();
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
{
double comp_moles = 0.0;
cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j];
if (comp_ref.Get_phase_name().size() == 0) continue;
/* First find exchange master species */
cxxNameDouble nd = comp_ref.Get_totals();
cxxNameDouble::iterator kit = nd.begin();
bool found_exchange = false;
for (; kit != nd.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr == NULL || elt_ptr->master == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
kit->first.c_str());
error_msg(error_string, CONTINUE);
continue;
}
if (elt_ptr->master->type == EX)
{
comp_moles = kit->second;
found_exchange = true;
}
}
//if (!found_exchange)
//{
// input_error++;
// error_string = sformatf(
// "Exchange formula does not contain an exchange master species, %s",
// comp_ref.Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n);
/* Now find the mineral on which exchanger depends... */
if (pp_assemblage_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Equilibrium_phases %d must be defined to use exchange related to mineral phase, %s",
n, comp_ref.Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
std::map<std::string, cxxPPassemblageComp>::iterator jit;
jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin();
for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++)
{
if (strcmp_nocase(comp_ref.Get_phase_name().c_str(), jit->first.c_str()) == 0)
{
break;
}
}
if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end())
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d",
comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
comp_ref.Set_phase_name(jit->first.c_str());
/* make exchanger concentration proportional to mineral ... */
conc = jit->second.Get_moles() * comp_ref.Get_phase_proportion();
if (found_exchange && comp_moles > 0.0)
{
/* parse formula */
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, 1.0);
free_check_null(temp_formula);
}
cxxNameDouble nd_formula = elt_list_NameDouble();
double comp_coef = 0;
for (kit = nd_formula.begin(); kit != nd_formula.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr->master->type == EX)
{
comp_coef = kit->second;
}
}
comp_ref.multiply(comp_coef * conc / comp_moles);
}
else /* comp_moles is zero, need to redefine totals from scratch */
{
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comp_ref.Set_totals(elt_list_NameDouble());
/*
* make sure exchange elements are in phase
*/
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, -comp_ref.Get_phase_proportion());
free_check_null(temp_formula);
}
int l;
struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE);
if (phase_ptr != NULL)
{
char* temp_formula = string_duplicate(phase_ptr->formula);
ptr = temp_formula;
get_elts_in_species(&ptr, 1.0);
free_check_null(temp_formula);
}
else
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d",
comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
qsort(elt_list, (size_t)count_elts,
(size_t)sizeof(struct elt_list), elt_list_compare);
elt_list_combine();
for (jj = 0; jj < count_elts; jj++)
{
if (elt_list[jj].elt->primary->s->type != EX
&& elt_list[jj].coef < 0)
{
input_error++;
error_string = sformatf(
"Stoichiometry of exchanger, %s * %g mol sites/mol phase,\n\tmust be a subset of the related phase %s, %s.",
comp_ref.Get_formula().c_str(),
(double)comp_ref.Get_phase_proportion(),
phase_ptr->name,
phase_ptr->formula);
error_msg(error_string, CONTINUE);
break;
}
}
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_min_surface(void)
/* ---------------------------------------------------------------------- */
/*
@ -3853,7 +4206,6 @@ tidy_min_surface(void)
assert(false);
}
cxxSurface *surface_ptr = &(kit->second);
if (!surface_ptr->Get_new_def()) continue;
if (!surface_ptr->Get_new_def())
continue;
if (surface_ptr->Get_n_user() < 0)
@ -4084,6 +4436,151 @@ tidy_min_surface(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_min_surface(void)
/* ---------------------------------------------------------------------- */
/*
* If surface is related to mineral, surface amount is
* set in proportion
*/
{
std::map<int, cxxSurface>::iterator kit;
for (kit = Rxn_surface_map.begin(); kit != Rxn_surface_map.end(); kit++)
{
cxxSurface* surface_ptr = &(kit->second);
if (surface_ptr->Get_n_user() < 0) continue;
for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++)
{
double comp_moles = 0.0;
cxxSurfaceComp* surface_comp_ptr = &(surface_ptr->Get_surface_comps()[j]);
if (surface_comp_ptr->Get_phase_name().size() == 0) continue;
cxxSurfaceCharge* surface_charge_ptr = NULL;
if (surface_ptr->Get_type() != cxxSurface::NO_EDL)
{
surface_charge_ptr = surface_ptr->Find_charge(surface_comp_ptr->Get_charge_name());
if (surface_charge_ptr == NULL)
{
input_error++;
error_string = sformatf("Data structure for surface charge not found "
"for %s ",
surface_comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
}
int n = surface_ptr->Get_n_user();
/* First find surface master species */
cxxNameDouble::iterator it;
for (it = surface_comp_ptr->Get_totals().begin(); it != surface_comp_ptr->Get_totals().end(); it++)
{
/* Find master species */
struct element* elt_ptr = element_store(it->first.c_str());
struct master* master_ptr = elt_ptr->master;
if (master_ptr == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
elt_ptr->name);
error_msg(error_string, CONTINUE);
continue;
}
if (master_ptr->type != SURF) continue;
comp_moles = it->second;
surface_comp_ptr->Set_master_element(elt_ptr->name);
break;
}
//if (surface_comp_ptr->Get_master_element().size() == 0)
//{
// input_error++;
// error_string = sformatf(
// "Surface formula does not contain a surface master species, %s",
// surface_comp_ptr->Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
/* Now find the mineral on which surface depends... */
cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n);
if (pp_assemblage_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Equilibrium_phases %d must be defined to use surface related to mineral phase, %s",
n, surface_comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
std::map<std::string, cxxPPassemblageComp>::iterator jit;
jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin();
for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++)
{
if (strcmp_nocase(surface_comp_ptr->Get_phase_name().c_str(),
jit->first.c_str()) == 0)
{
break;
}
}
if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end())
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to surface, %s, not found in Equilibrium_Phases %d",
surface_comp_ptr->Get_phase_name().c_str(), surface_comp_ptr->Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
int l;
struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE);
if (phase_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to surface, %s, not found in database.",
jit->first.c_str(), surface_comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
surface_comp_ptr->Set_phase_name(phase_ptr->name);
/* make surface concentration proportional to mineral ... */
LDBLE conc = jit->second.Get_moles() * surface_comp_ptr->Get_phase_proportion();
double grams = 0.0;
if (surface_charge_ptr != NULL)
{
grams = surface_charge_ptr->Get_grams();
}
if (comp_moles > 0.0)
{
surface_comp_ptr->multiply(conc / comp_moles);
}
else /* need to generate from scratch */
{
char* temp_formula = string_duplicate(surface_comp_ptr->Get_formula().c_str());
char* ptr = temp_formula;
count_elts = 0;
paren_count = 0;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
cxxNameDouble nd = elt_list_NameDouble();
surface_comp_ptr->Set_totals(nd);
}
if (grams > 0.0)
{
surface_charge_ptr->multiply(jit->second.Get_moles() / grams);
}
else if (surface_charge_ptr != NULL) /* need to generate from scratch */
{
surface_charge_ptr->Set_grams(jit->second.Get_moles());
surface_charge_ptr->Set_charge_balance(0.0);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_kin_surface(void)
/* ---------------------------------------------------------------------- */
/*
@ -4367,6 +4864,144 @@ tidy_kin_surface(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_kin_surface(void)
/* ---------------------------------------------------------------------- */
/*
* If surface is related to mineral, surface amount is
* set in proportion. Need to update surface if
* moles of kinetic reaction changes
*/
{
cxxKinetics* kinetics_ptr;
std::map<int, cxxSurface>::iterator it;
for (it = Rxn_surface_map.begin(); it != Rxn_surface_map.end(); it++)
{
cxxSurface* surface_ptr = &(it->second);
if (surface_ptr->Get_n_user() < 0) continue;
int n = surface_ptr->Get_n_user();
for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++)
{
double comp_moles = 0.0;
cxxSurfaceComp* comp_ptr = &(surface_ptr->Get_surface_comps()[j]);
if (comp_ptr->Get_rate_name().size() == 0) continue;
comp_ptr->Set_master_element("");
/* First find surface master species */
int k;
cxxNameDouble::iterator kit;
for (kit = comp_ptr->Get_totals().begin(); kit != comp_ptr->Get_totals().end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
struct master* master_ptr = elt_ptr->master;
if (master_ptr == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
elt_ptr->name);
error_msg(error_string, CONTINUE);
continue;
}
if (master_ptr->type != SURF) continue;
comp_ptr->Set_master_element(elt_ptr->name);
comp_moles = kit->second;
break;
}
if (comp_ptr->Get_master_element().size() == 0)
{
input_error++;
error_string = sformatf(
"Surface formula does not contain a surface master species, %s",
comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
/* Now find the kinetic reaction on which surface depends... */
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n)) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use surface related to kinetic reaction, %s",
n, comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
for (k = 0; k < (int)kinetics_ptr->Get_kinetics_comps().size(); k++)
{
cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
if (strcmp_nocase
(comp_ptr->Get_rate_name().c_str(),
kin_comp_ptr->Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == (int)kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to surface, %s, not found in Kinetics %d",
comp_ptr->Get_rate_name().c_str(), comp_ptr->Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
/* use database name for rate */
comp_ptr->Set_rate_name(kin_comp_ptr->Get_rate_name().c_str());
cxxSurfaceCharge* charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name());
if (surface_ptr->Get_type() != cxxSurface::NO_EDL)
{
charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name());
if (charge_ptr == NULL)
{
input_error++;
error_string = sformatf("Data structure for surface charge not found "
"for %s ",
comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
}
/* make surface concentration proportional to mineral ... */
LDBLE conc = kin_comp_ptr->Get_m() * comp_ptr->Get_phase_proportion();
double grams = 0.0;
if (charge_ptr != NULL) charge_ptr->Get_grams();
if (comp_moles > 0.0)
{
comp_ptr->multiply(conc / comp_moles);
}
else /* need to generate from scratch */
{
char* temp_formula = string_duplicate(comp_ptr->Get_formula().c_str());
char* ptr = temp_formula;
count_elts = 0;
paren_count = 0;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
cxxNameDouble nd = elt_list_NameDouble();
comp_ptr->Set_totals(nd);
}
if (grams > 0.0)
{
charge_ptr->multiply(kin_comp_ptr->Get_m() / grams);
}
else if (charge_ptr != NULL) /* need to generate from scratch */
{
charge_ptr->Set_grams(kin_comp_ptr->Get_m());
charge_ptr->Set_charge_balance(0.0);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
ss_prep(LDBLE t, cxxSS *ss_ptr, int print)
/* ---------------------------------------------------------------------- */
{

371
phreeqcpp/transport.cpp
View File

@ -1,4 +1,4 @@
#include "Utils.h"
#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "Exchange.h"
@ -585,7 +585,6 @@ transport(void)
if (overall_iterations > max_iter)
max_iter = overall_iterations;
cell_no = i;
mixrun = j;
if (multi_Dflag)
sprintf(token,
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@ -746,7 +745,6 @@ transport(void)
if (i == first_c && count_cells > 1)
kin_time /= 2;
cell_no = i;
mixrun = 0;
if (multi_Dflag)
sprintf(token,
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@ -794,6 +792,7 @@ transport(void)
j = 1;
for (; j <= nmix; j++) // loop on j
{
mixrun = j;
if (multi_Dflag && j == nmix && (transport_step % print_modulus == 0))
{
sprintf(token,
@ -850,7 +849,6 @@ transport(void)
if (overall_iterations > max_iter)
max_iter = overall_iterations;
cell_no = i;
mixrun = j;
if (multi_Dflag)
sprintf(token,
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
@ -2827,7 +2825,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
}
current_A = current_x / DDt * F_C_MOL;
for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++)
for (i = ifirst; i <= ilast + stagnant + ((bcon_last == 2 || (dV_dcell && stagnant)) ? 1 : 0); i++)
{
if (i <= ilast + 1)
{
@ -2879,10 +2877,12 @@ diffuse_implicit(LDBLE DDt, int stagnant)
for (icell = if1; icell != il1; icell += incr)
{
min_mol = min_dif_M * ct[icell].kgw;
if (min_mol < 1e-13)
min_mol = 1e-13;
dum1 = dum2 = 0;
sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1);
if (stagnant)
if (stagnant && mixf_stag[icell][cp])
{
i1 = (icell == 0 ? c1 + 1 : icell == c1 ? cc1 : icell + c1);
sptr_stag = Utilities::Rxn_find(Rxn_solution_map, i1);
@ -2890,6 +2890,13 @@ diffuse_implicit(LDBLE DDt, int stagnant)
else
sptr_stag = NULL;
//if (!cp)
//{
// ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0;
// if (sptr_stag)
// ct[i1].J_ij_sum = 0.0;
//}
if (!strcmp(ct[icell].m_s[cp].name, "H"))
{
dummy = ct[icell].m_s[cp].tot1;
@ -2933,7 +2940,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
}
sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy);
}
//if (cp == count_m_s - 1) // transport the charge imbalance (but doesn't improve the change of H2 or O2).
//if (cp == count_m_s - 1) // transport the charge imbalance
//{
// sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
// sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
@ -2942,180 +2949,166 @@ diffuse_implicit(LDBLE DDt, int stagnant)
//}
continue;
}
// see if icell - incr has negative moles, subtract it, so that it is canceled...
if (icell > 0 && icell <= ilast &&
(it1 = neg_moles.find(icell - incr)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
{
ct[icell].m_s[cp].tot1 += it2->second;
neg_moles.erase(it1);
els.erase(it2);
neg_moles.insert(std::make_pair(icell - incr, els));
}
dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name];
if (!dum1)
dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name];
if (sptr_stag)
dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name];
// 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)))
{
dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name);
if (dum1)
sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1;
}
dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name];
if (!dum2)
{
dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
if (dum2)
sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
}
if (sptr_stag)
{
dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name];
if (!dum_stag)
{
dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name);
if (dum_stag)
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag;
}
}
// check for negative moles, add moles from the donnan layer when necessary and available...
if (ct[icell].m_s[cp].tot1 > dum1 &&
(dV_dcell || (icell > 0 && icell <= ilast)))
{
if (!ct[icell].dl_s)
ct[icell].m_s[cp].tot1 = dum1;
else
if (ct[icell].dl_s > 1e-8)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell);
if (s_ptr)
{
dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 - dum1 + min_mol);
dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 + ct[icell].m_s[cp].tot_stag - dum1 + min_mol);
}
}
sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1;
if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol)
ct[icell].m_s[cp].tot1 = dum1 - ct[icell].m_s[cp].tot_stag - min_mol;
}
// check for negative moles, in the other cell...
ct[icell].m_s[cp].tot2 = ct[icell].m_s[cp].tot1;
if (-ct[icell].m_s[cp].tot2 > dum2 &&
(dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
dum = 0;
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)*/))
{
if (!ct[icell + 1].dl_s)
ct[icell].m_s[cp].tot2 = -dum2;
else
dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
if (ct[icell + 1].dl_s > 1e-8)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1);
if (s_ptr)
{
dum2 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(ct[icell].m_s[cp].tot2 + dum2 - min_mol));
}
sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
if (-ct[icell].m_s[cp].tot2 > dum2)
ct[icell].m_s[cp].tot2 = -dum2;
}
sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol)
ct[icell].m_s[cp].tot2 = -dum2 + min_mol + dum;
}
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 (!cp)
{
ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0;
if (sptr_stag)
ct[i1].J_ij_sum = 0.0;
}
if (dV_dcell || (icell > 0 && icell <= ilast))
{
dum1 -= ct[icell].m_s[cp].tot1;
if (sptr_stag)
dum1 -= ct[icell].m_s[cp].tot_stag;
sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 0e-16);
if (dum1 < -min_mol && incr > 0)
dum = ct[icell].m_s[cp].tot1;
if (stagnant)
dum += ct[icell].m_s[cp].tot_stag;
dum1 -= dum;
sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : min_mol);
if (dum1 < 0)
{
dum += dum1 - min_mol;
if ((it1 = neg_moles.find(icell)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum1 += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1));
neg_moles.erase(icell);
neg_moles.insert(std::make_pair(icell, els));
}
else
ct[icell].J_ij_sum += dum1 * ct[icell].m_s[cp].charge;
//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 || (icell >= 0 && icell < ilast)/* || (icell == ilast && bcon_last == 2)*/)
{
dum2 += ct[icell].m_s[cp].tot1;
sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 0e-16);
if (dum2 < -min_mol && incr < 0)
dum = ct[icell].m_s[cp].tot1;
if (stagnant && icell == c && dV_dcell)
dum -= ct[c1].m_s[cp].tot_stag;
dum2 += dum;
sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol);
if (dum2 < 0)
{
dum -= dum2 - min_mol;
if ((it1 = neg_moles.find(icell + 1)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum2 += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2));
neg_moles.erase(icell + 1);
neg_moles.insert(std::make_pair(icell + 1, els));
}
else
ct[icell + 1].J_ij_sum += dum2 * ct[icell].m_s[cp].charge;
//ct[icell + 1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
}
if (sptr_stag)
{
if ((it1 = neg_moles.find(i1)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
{
ct[icell].m_s[cp].tot_stag += it2->second;
// the negative moles are canceled...
neg_moles.erase(it1);
els.erase(it2);
neg_moles.insert(std::make_pair(i1, els));
}
dummy = ct[icell].m_s[cp].tot_stag;
dum = ct[icell].m_s[cp].tot_stag;
if (icell == c)
dum += ct[c1].m_s[cp].tot_stag;
if (dum_stag + dum < 0)
{
if (dV_dcell)
{
if ((dum = sptr2->Get_totals()[ct[icell + 1].m_s[cp].name] - ct[icell + 1].m_s[cp].tot_stag) > 0)
{
sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum;
dummy += ct[icell + 1].m_s[cp].tot_stag;
}
else
{
dummy += sptr2->Get_totals()[ct[icell].m_s[cp].name];
sptr2->Get_totals()[ct[icell].m_s[cp].name] = 0.0;
if (dum < -min_mol)
{
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum));
neg_moles.erase(icell);
neg_moles.insert(std::make_pair(icell, els));
}
}
}
else
dummy += ct[icell + 1].m_s[cp].tot_stag;
}
if (-dummy > dum_stag)
{
dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name);
if (ct[i1].dl_s)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, i1);
if (s_ptr)
{
dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dummy + dum_stag - min_mol));
dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dum + dum_stag - min_mol));
}
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag;
}
}
dum_stag += dummy;
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : 0e-16);
if (dum_stag < -min_mol)
dum_stag += dum;
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : min_mol);
if (dum_stag < 0)
{
dum -= dum_stag - min_mol;
if ((it1 = neg_moles.find(i1)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum_stag += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum_stag));
neg_moles.erase(i1);
neg_moles.insert(std::make_pair(i1, els));
}
else
ct[i1].J_ij_sum += dum_stag * ct[icell].m_s[cp].charge;
//ct[i1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
}
//if (cp == count_m_s - 1)
// reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted...
if (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)
{
sptr1 = Utilities::Rxn_find(Rxn_solution_map, 1);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, 2);
dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[1].kgw * (1 - ct[1].dl_s) - dummy;
dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[2].kgw * (1 - ct[2].dl_s) - dummy;
if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
{
dum = cell_data[1].mid_cell_x / cell_data[2].mid_cell_x;
//ct[1].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[1].m_s[cp].charge;
dum1 = (dummy + dum * dum2) * ct[1].kgw / (1 - ct[1].dl_s);
sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
//ct[1].J_ij_sum += dum1 * ct[1].m_s[cp].charge;
}
}
dummy = Utilities::Rxn_find(Rxn_solution_map, c1)->Get_totals()[ct[0].m_s[cp].name] / ct[c1].kgw * (1 - ct[c1].dl_s);
if (dummy > 1e-6)
{
sptr1 = Utilities::Rxn_find(Rxn_solution_map, c);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, c_1);
dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[c].kgw * (1 - ct[c].dl_s) - dummy;
dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[c_1].kgw * (1 - ct[c_1].dl_s) - dummy;
if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
{
dum = (cell_data[c].mid_cell_x - cell_data[c_1].mid_cell_x) /
(cell_data[c1].mid_cell_x - cell_data[c_1].mid_cell_x);
//ct[c].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[c].m_s[cp].charge;
dum1 = (dummy + (1 - dum) * dum2) * ct[c].kgw / (1 - ct[c].dl_s);
sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
//ct[c].J_ij_sum += dum1 * ct[c].m_s[cp].charge;
}
}
}
//if (cp == count_m_s - 1)
//{
// sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
// sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
@ -3664,6 +3657,46 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant)
}
return (OK);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant)
/* ---------------------------------------------------------------------- */
{
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j));
// At filterends, concentrations of ions change step-wise to the DL.
// We take the harmonic mean for f_free, the average for the DL.
if (ct[icell].v_m[k].z)
{
if (!g_i && g_j)
{
ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) +
b_i * (1 - free_j) / 4 + b_j * g_j / 4;
}
else if (g_i && !g_j)
ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) +
b_j * (1 - free_i) / 4 + b_i * g_i / 4;
}
// for boundary cells...
if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2;
}
else
{
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j);
else if (icell == count_cells && jcell == count_cells + 1)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i);
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
return;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
@ -3848,7 +3881,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
if (dl_aq1 > 0)
ct[icell].dl_s = dl_aq1 / t_aq1;
if (dl_aq2 > 0)
ct[icell].dl_s = ct[jcell].dl_s = dl_aq2 / t_aq2;
{
ct[jcell].dl_s = dl_aq2 / t_aq2;
if (!ct[icell].dl_s)
ct[icell].dl_s = 1e-8; // used in implicit appt
}
if (il_calcs)
{
@ -4142,7 +4179,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else
{
dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
dum1 = it_sc->Get_mass_water() / t_aq2;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_j += pow(dum2, ct[icell].v_m[k].z) * dum1;
}
@ -4152,32 +4189,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
}
b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * (f_free_j + g_j / ct[icell].visc2);
if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
else if (icell == all_cells - 1 && stagnant)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf *= 2 for this 'reservoir' cell in the input */
b_i = A1 * sol_D[icell].spec[i].Dwt;
b_j = A2;
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_j *= sol_D[icell].spec[i].Dwt;
else
{
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_j *= sol_D[icell].spec[i].Dwt;
else
{
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
dum2 *= sol_D[jcell].viscos_f;
b_j *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == 0 && !stagnant)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == 3 && stagnant && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for stagnant cell 3 in the input */
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
dum2 *= sol_D[jcell].viscos_f;
b_j *= dum2;
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
}
@ -4249,7 +4272,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else
{
dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
dum1 = it_sc->Get_mass_water() / t_aq1;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_i += pow(dum2, ct[icell].v_m[k].z) * dum1;
}
@ -4266,31 +4289,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j *= sol_D[jcell].spec[j].erm_ddl;
}
}
b_i = A1 * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
if (icell == 0 && !stagnant)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == 3 && stagnant && g_j && !g_i)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for 'reservoir' cell 3 in the input */
b_i = A1;
b_j = A2 * sol_D[jcell].spec[j].Dwt;
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_i *= sol_D[jcell].spec[j].Dwt;
else
{
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_i *= sol_D[jcell].spec[j].Dwt;
else
{
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
dum2 *= sol_D[icell].viscos_f;
b_i *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
else if (jcell == all_cells - 1 && stagnant && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf * 2 for this 'reservoir' cell in the input */
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
dum2 *= sol_D[icell].viscos_f;
b_i *= dum2;
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
}
@ -4373,28 +4383,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j *= sol_D[jcell].spec[j].erm_ddl;
}
}
b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
// but for boundary cells...
if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2;
}
else
{
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
ct[icell].v_m[k].b_ij = b_j;
else if (icell == count_cells && jcell == count_cells + 1)
ct[icell].v_m[k].b_ij = b_i;
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
b_i = A1 * sol_D[icell].spec[i].Dwt;
b_j = A2 * sol_D[jcell].spec[j].Dwt;
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
//ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; // appt: this could give an incorrect large factor for implicit
//if (fabs(ddlm) > 1e-10)