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 '710d0096aa62355bf6c212efb3f6f9392ccd3f37'

Browse Source
This commit is contained in:
Darth Vader 2021-03-12 16:57:53 +00:00
commit 3d05ce1743
12 changed files with 262 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
phreeqcpp/PBasic.cpp
View File

@ -1387,6 +1387,22 @@ listtokens(FILE * f, tokenrec * l_buf)
output_msg("LK_SPECIES");
break;
case tokdelta_h_species:
output_msg("DELTA_H_SPECIES");
break;
case tokdelta_h_phase:
output_msg("DELTA_H_PHASE");
break;
case tokdh_a0:
output_msg("DH_A0");
break;
case tokdh_bdot:
output_msg("DH_BDOT");
break;
case toklk_named:
output_msg("LK_NAMED");
break;
@ -2457,6 +2473,31 @@ factor(struct LOC_exec * LINK)
}
break;
case tokdelta_h_phase:
{
const char* str = stringfactor(STR1, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_deltah_p(str);
}
break;
case tokdelta_h_species:
{
const char* str = stringfactor(STR1, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_deltah_s(str);
}
break;
case tokdh_a0:
{
const char* str = stringfactor(STR1, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->dh_a0(str);
}
break; case tokdh_bdot:
{
const char* str = stringfactor(STR1, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->dh_bdot(str);
}
break;
case toklk_named:
{
const char* str = stringfactor(STR1, LINK);
@ -7416,7 +7457,11 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("lk_species", PBasic::toklk_species),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("lk_named", PBasic::toklk_named),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("lk_phase", PBasic::toklk_phase),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("delta_h_phase", PBasic::tokdelta_h_phase),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("delta_h_species", PBasic::tokdelta_h_species),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sum_species", PBasic::toksum_species),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_a0", PBasic::tokdh_a0),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("dh_bdot", PBasic::tokdh_bdot),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sum_gas", PBasic::toksum_gas),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sum_s_s", PBasic::toksum_s_s),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("calc_value", PBasic::tokcalc_value),

4
phreeqcpp/PBasic.h
View File

@ -259,6 +259,10 @@ public:
toklk_species,
toklk_named,
toklk_phase,
tokdelta_h_phase,
tokdelta_h_species,
tokdh_a0,
tokdh_bdot,
toksum_species,
toksum_gas,
toksum_s_s,

6
phreeqcpp/Phreeqc.h
View File

@ -114,6 +114,10 @@ public:
LDBLE calc_logk_n(const char *name);
LDBLE calc_logk_p(const char *name);
LDBLE calc_logk_s(const char *name);
LDBLE calc_deltah_s(const char* name);
LDBLE calc_deltah_p(const char* name);
LDBLE dh_a0(const char* name);
LDBLE dh_bdot(const char* name);
LDBLE calc_surface_charge(const char *surface_name);
LDBLE calc_t_sc(const char *name);
LDBLE diff_layer_total(const char *total_name, const char *surface_name);
@ -1701,7 +1705,7 @@ protected:
bool output_newline;
inline void Set_output_newline(bool tf) { this->output_newline = tf;}
inline bool Get_output_newline() { return this->output_newline;}
LDBLE *llnl_temp, *llnl_adh, *llnl_bdh, *llnl_bdot, *llnl_co2_coefs, a_llnl, b_llnl;
LDBLE *llnl_temp, *llnl_adh, *llnl_bdh, *llnl_bdot, *llnl_co2_coefs, a_llnl, b_llnl, bdot_llnl;
int llnl_count_temp, llnl_count_adh, llnl_count_bdh, llnl_count_bdot,
llnl_count_co2_coefs;

1
phreeqcpp/SelectedOutput.h
View File

@ -181,7 +181,6 @@ protected:
bool charge_balance;
bool percent_error;
bool new_line;
//bool punch_newline;
// as-is set flags
bool set_user_punch;

2
phreeqcpp/SurfaceComp.cxx
View File

@ -432,7 +432,7 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
{
std::ostringstream oss;
oss <<
"Cannot mix surface components related to phase with exchange components related to kinetics, "
"Cannot mix surface components related to phase with surface components related to kinetics, "
<< this->formula;
error_msg(oss.str().c_str(), CONTINUE);
return;

112
phreeqcpp/basicsubs.cpp
View File

@ -729,6 +729,117 @@ calc_logk_s(const char *name)
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
dh_a0(const char* name)
/* ---------------------------------------------------------------------- */
{
char token[MAX_LENGTH];
struct species* s_ptr;
double a = -999.99;
strcpy(token, name);
s_ptr = s_search(token);
if (s_ptr != NULL)
{
a = s_ptr->dha;
}
return (a);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
dh_bdot(const char* name)
/* ---------------------------------------------------------------------- */
{
char token[MAX_LENGTH];
struct species* s_ptr;
double b = -999.99;
if (llnl_count_temp > 0)
{
b = bdot_llnl;
}
else
{
strcpy(token, name);
s_ptr = s_search(token);
if (s_ptr != NULL)
{
b = s_ptr->dhb;
}
}
return (b);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_deltah_p(const char* name)
/* ---------------------------------------------------------------------- */
{
int i, j;
char token[MAX_LENGTH];
struct phase* phase_ptr;
LDBLE lkm, lkp;
LDBLE l_logk[MAX_LOG_K_INDICES];
double dh = -999.99;
strcpy(token, name);
phase_ptr = phase_bsearch(token, &j, FALSE);
if (phase_ptr != NULL)
{
struct reaction* reaction_ptr;
if (phase_ptr->replaced)
reaction_ptr = phase_ptr->rxn_s;
else
reaction_ptr = phase_ptr->rxn;
/*
* Print saturation index
*/
reaction_ptr->logk[delta_v] = calc_delta_v(reaction_ptr, true) -
phase_ptr->logk[vm0];
if (reaction_ptr->logk[delta_v])
mu_terms_in_logk = true;
for (i = 0; i < MAX_LOG_K_INDICES; i++)
{
l_logk[i] = 0.0;
}
//lk = k_calc(reaction_ptr->logk, tk_x, patm_x * PASCAL_PER_ATM);
select_log_k_expression(reaction_ptr->logk, l_logk);
add_other_logk(l_logk, phase_ptr->count_add_logk, phase_ptr->add_logk);
lkm = k_calc(l_logk, tk_x - 1.0, patm_x * PASCAL_PER_ATM);
lkp = k_calc(l_logk, tk_x + 1.0, patm_x * PASCAL_PER_ATM);
dh = (lkp - lkm) / 2.0 * LOG_10 * R_KJ_DEG_MOL * pow(tk_x, 2.0);
}
return (dh);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_deltah_s(const char* name)
/* ---------------------------------------------------------------------- */
{
int i;
char token[MAX_LENGTH];
struct species* s_ptr;
LDBLE lkm, lkp, l_logk[MAX_LOG_K_INDICES];
double dh = -999.99;
strcpy(token, name);
s_ptr = s_search(token);
if (s_ptr != NULL)
{
/* calculate delta_v for the reaction... */
s_ptr->logk[delta_v] = calc_delta_v(s_ptr->rxn, false);
for (i = 0; i < MAX_LOG_K_INDICES; i++)
{
l_logk[i] = 0.0;
}
select_log_k_expression(s_ptr->logk, l_logk);
mu_terms_in_logk = true;
add_other_logk(l_logk, s_ptr->count_add_logk, s_ptr->add_logk);
lkm = k_calc(l_logk, tk_x-1.0, patm_x * PASCAL_PER_ATM);
lkp = k_calc(l_logk, tk_x + 1.0, patm_x * PASCAL_PER_ATM);
dh = (lkp - lkm) / 2.0 * LOG_10 * R_KJ_DEG_MOL * pow(tk_x,2.0);
return (dh);
}
return (0.0);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_surface_charge(const char *surface_name)
/* ---------------------------------------------------------------------- */
{
@ -738,7 +849,6 @@ calc_surface_charge(const char *surface_name)
LDBLE charge;
struct rxn_token_temp *token_ptr;
struct master *master_ptr;
/*
* Go through species, sum charge
*/

12
phreeqcpp/class_main.cpp
View File

@ -277,13 +277,13 @@ write_banner(void)
char buffer[80];
int len, indent;
screen_msg(
" █▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█\n");
" █▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█\n");
screen_msg(
" ║ ║\n");
/* version */
#ifdef NPP
len = sprintf(buffer, "* PHREEQC-%s *", "3.5.2");
len = sprintf(buffer, "* PHREEQC-%s *", "3.6.5");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -307,7 +307,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "August 1, 2019");
len = sprintf(buffer, "%s", "February 24, 2021");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -316,7 +316,7 @@ write_banner(void)
44 - indent - len, ' '));
screen_msg(
" █▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄█\n\n");
" █▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄█\n\n");
return 0;
}
@ -491,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));
@ -499,7 +499,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
output_msg(sformatf(" Input file: %s\n", in_file));
output_msg(sformatf(" Output file: %s\n", out_file));
#ifdef NPP
output_msg(sformatf("Using PHREEQC: version 3.5.2, compiled August 1, 2019\n"));
output_msg(sformatf("Using PHREEQC: version 3.6.5, compiled February 24, 2021\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token));
#ifdef NPP

80
phreeqcpp/inverse.cpp
View File

@ -2029,14 +2029,31 @@ print_model(struct inverse *inv_ptr)
(double) d3));
}
output_msg(sformatf( "\n%-25.25s %2s %12.12s %12.12s\n",
"Phase mole transfers:", " ", "Minimum", "Maximum"));
// appt, calculate and print SI's
int i1, i2;
LDBLE t_i, p_i, iap, lk, t;
const char *name;
struct rxn_token *rxn_ptr;
struct reaction *reaction_ptr;
output_msg(sformatf( "\n%-25.25s %2s %12.12s %12.12s %-18.18s (Approximate SI in solution ",
"Phase mole transfers:", " ", "Minimum", "Maximum", "Formula"));
for (i = 0; i < inv_ptr->count_solns - 1; i++)
output_msg(sformatf("%d, ", inv_ptr->solns[i]));
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
t_i = solution_ptr->Get_tc() + 273.15;
p_i = solution_ptr->Get_patm();
output_msg(sformatf("%d at %3d K, %3d atm)\n", inv_ptr->solns[i], int(t_i), int(floor(p_i + 0.5))));
p_i *= PASCAL_PER_ATM;
for (i = col_phases; i < col_redox; i++)
{
if (equal(inv_delta1[i], 0.0, toler) == TRUE &&
equal(min_delta[i], 0.0, toler) == TRUE &&
equal(max_delta[i], 0.0, toler) == TRUE)
continue;
d1 = inv_delta1[i];
d2 = min_delta[i];
d3 = max_delta[i];
@ -2047,11 +2064,62 @@ print_model(struct inverse *inv_ptr)
if (equal(d3, 0.0, MIN_TOTAL_INVERSE) == TRUE)
d3 = 0.0;
output_msg(sformatf(
"%15.15s %12.3e %12.3e %12.3e %s\n", col_name[i],
(double) d1, (double) d2, (double) d3,
inv_ptr->phases[i - col_phases].phase->formula));
}
"%15.15s %12.3e %12.3e %12.3e %-25.25s (", col_name[i],
(double)d1, (double)d2, (double)d3, inv_ptr->phases[i - col_phases].phase->formula));
i1 = 0;
for (; i1 < count_phases; i1++)
{
if (Utilities::strcmp_nocase(phases[i1]->name, col_name[i]))
continue;
reaction_ptr = phases[i1]->rxn_s;
for (i2 = 0; i2 < inv_ptr->count_solns; i2++)
{
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i2]);
reaction_ptr->logk[delta_v] = calc_delta_v(reaction_ptr, true) - phases[i1]->logk[vm0];
if (reaction_ptr->logk[delta_v])
mu_terms_in_logk = true;
lk = k_calc(reaction_ptr->logk, t_i, p_i);
iap = 0.0;
for (rxn_ptr = reaction_ptr->token + 1; rxn_ptr->s != NULL; rxn_ptr++)
{
t = 0;
if (rxn_ptr->s == s_eminus)
t = -solution_ptr->Get_pe();
else if (!Utilities::strcmp_nocase(rxn_ptr->s->name, "H2O"))
t = log10(solution_ptr->Get_ah2o());
else if (!Utilities::strcmp_nocase(rxn_ptr->s->name, "H+"))
t = -solution_ptr->Get_ph();
else
{
if (rxn_ptr->s->secondary)
name = rxn_ptr->s->secondary->elt->name;
else
name = rxn_ptr->s->primary->elt->name;
t = solution_ptr->Get_master_activity()[name];
}
if (t)
iap += t * rxn_ptr->coef;
else
{
iap = -999; break;
}
}
if (iap == -999)
output_msg(sformatf(" "));
else
output_msg(sformatf("%6.2f", iap - lk));
if (i2 < inv_ptr->count_solns - 1)
output_msg(sformatf(","));
}
}
output_msg(sformatf(")\n"));
}
output_msg(sformatf( "\n%-25.25s\n", "Redox mole transfers:"));
for (i = col_redox; i < col_epsilon; i++)
{

8
phreeqcpp/kinetics.cpp
View File

@ -2005,7 +2005,13 @@ set_reaction(int i, int use_mix, int use_kinetics)
/*
* Find surface
*/
dl_type_x = cxxSurface::NO_DL;
if (use.Get_surface_in() && use.Get_kinetics_in() && use.Get_kinetics_ptr() && !use.Get_kinetics_ptr()->Get_use_cvode() && reaction_step > 1)
{
// use.Set_surface_ptr(Utilities::Rxn_find(Rxn_surface_map, i));
// appt: we may come here with zero kinetic reaction, but surface may have to keep DONNAN_DL
}
else
dl_type_x = cxxSurface::NO_DL;
if (use.Get_surface_in() == TRUE)
{
use.Set_surface_ptr(Utilities::Rxn_find(Rxn_surface_map, i));

2
phreeqcpp/model.cpp
View File

@ -568,7 +568,7 @@ gammas(LDBLE mu)
*/
int i, j;
int ifirst, ilast;
LDBLE f, bdot_llnl, log_g_co2, dln_g_co2, c2_llnl;
LDBLE f, log_g_co2, dln_g_co2, c2_llnl;
LDBLE c1, c2, a, b;
LDBLE muhalf, equiv;

2
phreeqcpp/prep.cpp
View File

@ -57,7 +57,7 @@ prep(void)
//if (!same_model && !switch_numerical)
// numerical_fixed_volume = false;
if (same_model == FALSE /*|| switch_numerical*/)
if (same_model == FALSE || !my_array/*|| switch_numerical*/)
{
clear();
setup_unknowns();

14
phreeqcpp/transport.cpp
View File

@ -1,4 +1,4 @@
#include "Utils.h"
#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "Exchange.h"
@ -2650,7 +2650,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
// decompose A in LU : store L in A[..][0..1] and U in A[..][2] ...
for (i = 1; i <= count_cells + 1; i++)
{
A[i - 1][2] = A[i - 1][2] / A[i - 1][1];
A[i - 1][2] /= A[i - 1][1];
A[i][1] -= A[i][0] * A[i - 1][2];
}
// solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1, Ct1 was put in Ct2 ...
@ -2786,7 +2786,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
}
dVc = j_0e * current_cells[ifirst].R;
cell_data[ifirst + 1].potV = cell_data[ifirst].potV + dVc;
for (i = ifirst + 1; i < ilast; i++)
for (i = ifirst + 1; i <= ilast; i++)
{
dVc = current_cells[i].R * (current_x - current_cells[i].dif);
//if (((dV_dcell && (dVc * j_0e > 0)) ||
@ -2825,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 || (dV_dcell && stagnant)) ? 1 : 0); i++)
for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++)
{
if (i <= ilast + 1)
{
@ -2977,7 +2977,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
if (icell == c && sptr_stag && ct[c1].m_s[cp].tot_stag)
dum = ct[c1].m_s[cp].tot_stag;
if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol &&
(dV_dcell || (icell >= 0 && icell <= ilast)/* || (icell == ilast && bcon_last == 2)*/))
(dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
{
dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
if (ct[icell + 1].dl_s > 1e-8)
@ -3016,7 +3016,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
//ct[icell].J_ij_sum -= dum * ct[icell].m_s[cp].charge;
}
if (dV_dcell || (icell >= 0 && icell < ilast)/* || (icell == ilast && bcon_last == 2)*/)
if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
{
dum = ct[icell].m_s[cp].tot1;
if (stagnant && icell == c && dV_dcell)
@ -3072,7 +3072,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
}
// reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted...
if (icell == il1 - incr && ct[0].m_s != NULL && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1)
if (dV_dcell && 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)