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 'f90fc5afeaa5347ef057416551bb406e8ce5d855'

Browse Source
This commit is contained in:
Scott R Charlton 2019-01-28 18:20:25 -07:00
commit f56958652e
11 changed files with 171 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
phreeqcpp/PBasic.cpp
View File

@ -1777,7 +1777,11 @@ void PBasic::
_ASSERTE(nIDErrPrompt == 0);
nIDErrPrompt = IDS_ERR_SYNTAX;
}
errormsg(strcat(str, l_s));
strcat(str, l_s);
strcat(str, " in line: ");
if (strcmp(inbuf, "run"))
strcat(str, inbuf);
errormsg(str);
}
void PBasic::
@ -1790,7 +1794,11 @@ void PBasic::
_ASSERTE(nIDErrPrompt == 0);
nIDErrPrompt = IDS_ERR_MISMATCH;
}
errormsg(strcat(str, l_s));
strcat(str, l_s);
strcat(str, " in line: ");
if (strcmp(inbuf, "run"))
strcat(str, inbuf);
errormsg(str);
}
void PBasic::
@ -6031,7 +6039,10 @@ exec(void)
_ASSERTE(nIDErrPrompt == 0);
nIDErrPrompt = IDS_ERR_ILLEGAL;
}
errormsg("Illegal command");
strcat(STR1, "Illegal command in line: ");
if (strcmp(inbuf, "run"))
strcat(STR1, inbuf);
errormsg(STR1);
break;
}
}

4
phreeqcpp/Phreeqc.h
View File

@ -493,6 +493,7 @@ public:
int calc_fixed_volume_gas_pressures(void);
int calc_ss_fractions(void);
int gammas(LDBLE mu);
int gammas_a_f(int i);
int initial_guesses(void);
int revise_guesses(void);
int ss_binary(cxxSS *ss_ptr);
@ -535,7 +536,8 @@ public:
struct theta_param *theta_param_alloc(void);
int theta_param_init(struct theta_param *theta_param_ptr);
void pitzer_make_lists(void);
int gammas_pz(void);
//int gammas_pz(void);
int gammas_pz(bool exch_a_f);
int model_pz(void);
int pitzer(void);
int pitzer_clean_up(void);

10
phreeqcpp/class_main.cpp
View File

@ -277,8 +277,8 @@ write_banner(void)
/* version */
#ifdef NPP
//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.6 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.4.6");
//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -302,7 +302,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "June 20, 2018");
len = sprintf(buffer, "%s", "January 17, 2019");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -492,12 +492,12 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
screen_msg(sformatf("Database file: %s\n\n", token));
strcpy(db_file, token);
#ifdef NPP
//output_msg(sformatf("Using PHREEQC: version 3.4.6 Ampère, compiled June 20, 2018\n"));
//output_msg(sformatf("Using PHREEQC: version 3.4.8 Ampère, compiled January 17, 2019\n"));
#endif
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.4.6, compiled June 20, 2018\n"));
output_msg(sformatf("Using PHREEQC: version 3.4.8, compiled January 17, 2019\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token));
#ifdef NPP

15
phreeqcpp/kinetics.cpp
View File

@ -2474,6 +2474,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
*/
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, i);
if (kin_time <= 0 ||
(kinetics_ptr && kinetics_ptr->Get_kinetics_comps().size() == 0) ||
(state == REACTION && use.Get_kinetics_in() == FALSE) ||
(state == TRANSPORT && kinetics_ptr == NULL) ||
(state == PHAST && kinetics_ptr == NULL) ||
@ -2661,10 +2662,13 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
while (flag != SUCCESS)
{
sum_t += cvode_last_good_time;
error_string = sformatf(
if (state != TRANSPORT)
{
error_string = sformatf(
"CVode incomplete at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d, continuing...\n",
(int) iopt[NST], cell_no, (double) sum_t, m_iter + 1);
warning_msg(error_string);
(int)iopt[NST], cell_no, (double)sum_t, m_iter + 1);
warning_msg(error_string);
}
#ifdef DEBUG_KINETICS
if (m_iter > 5)
dump_kinetics_stderr(cell_no);
@ -2675,7 +2679,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
{
m_temp = (LDBLE *) free_check_null(m_temp);
m_original = (LDBLE *) free_check_null(m_original);
error_msg("Repeated restart of integration.", STOP);
error_string = sformatf(
"Restart of integration at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d.",
(int)iopt[NST], cell_no, (double)sum_t, m_iter - 1);
error_msg(error_string, STOP);
}
tout1 = tout - sum_t;
t = 0;

160
phreeqcpp/model.cpp
View File

@ -563,14 +563,14 @@ gammas(LDBLE mu)
*/
int i, j;
int ifirst, ilast;
LDBLE d1, d2, d3, f, a_llnl, b_llnl, bdot_llnl, log_g_co2, dln_g_co2, c2_llnl;
LDBLE f, a_llnl, b_llnl, bdot_llnl, log_g_co2, dln_g_co2, c2_llnl;
LDBLE c1, c2, a, b;
LDBLE muhalf, equiv;
/* Initialize */
if (mu <= 0) mu = 1e-10;
if (pitzer_model == TRUE)
return gammas_pz();
return gammas_pz(true);
if (sit_model == TRUE)
return gammas_sit();
a_llnl = b_llnl = bdot_llnl = log_g_co2 = dln_g_co2 = c2_llnl = 0;
@ -683,57 +683,58 @@ gammas(LDBLE mu)
/*
* Find CEC
* z contains charge of cation for exchange species, alk contains cec
* correct activity for Gapon-type exchange eqns: Ca0.5X uses (gamma_Ca)^0.5
*/
/* !!!!! */
if (calculating_deriv)
continue;
LDBLE coef, z;
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
if (s_x[i]->rxn_x->token[j].s->type == EX)
{
s_x[i]->alk =
s_x[i]->rxn_x->token[j].s->primary->unknown->moles;
break;
//break;
}
else if (s_x[i]->rxn_x->token[j].s->type <= HPLUS)
{
coef = s_x[i]->rxn_x->token[j].coef;
z = s_x[i]->rxn_x->token[j].s->z;
}
}
if (s_x[i]->exch_gflag == 1 && s_x[i]->alk > 0)
if (!use.Get_exchange_ptr()->Get_pitzer_exchange_gammas())
{
if (s_x[i]->primary != NULL)
{
s_x[i]->lg = 0.0;
s_x[i]->dg = 0.0;
}
else
{
if (s_x[i]->alk <= 0)
s_x[i]->lg = 0.0;
else
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = 0.0;
}
}
else if (s_x[i]->exch_gflag == 1 && s_x[i]->alk > 0)
{
/* Davies */
d1 = s_x[i]->lg;
s_x[i]->lg = -s_x[i]->equiv * s_x[i]->equiv * a *
s_x[i]->lg = -coef * z * z * a *
(muhalf / (1.0 + muhalf) - 0.3 * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
if (s_x[i]->a_f && s_x[i]->primary == NULL)
{
d2 = s_x[i]->moles * s_x[i]->equiv / s_x[i]->alk;
if (d2 > 1) d2 = 1;
d2 = s_x[i]->lg - s_x[i]->a_f * (1 - d2);
d3 = 0.89;
if (iterations < 10) d3 = 0.7; else d3 = 0.89;
s_x[i]->lg = d3 * d1 + (1 - d3) * d2;
}
s_x[i]->dg =
c1 * s_x[i]->equiv * s_x[i]->equiv * s_x[i]->moles;
c1 * coef * z * z * s_x[i]->moles;
}
else if (s_x[i]->exch_gflag == 2 && s_x[i]->alk > 0)
{
/* Extended D-H, WATEQ D-H */
d1 = s_x[i]->lg;
s_x[i]->lg = -a * muhalf * s_x[i]->equiv * s_x[i]->equiv /
(1.0 + s_x[i]->dha * b * muhalf) + s_x[i]->dhb * mu +
s_x[i]->lg = coef * (-a * muhalf * z * z /
(1.0 + s_x[i]->dha * b * muhalf) + s_x[i]->dhb * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
if (s_x[i]->a_f && s_x[i]->primary == NULL)
{
d2 = s_x[i]->moles * s_x[i]->equiv / s_x[i]->alk;
if (d2 > 1) d2 = 1;
d2 = s_x[i]->lg - s_x[i]->a_f * (1 - d2);
d3 = 0.89;
if (iterations < 10) d3 = 0.7; else d3 = 0.89;
s_x[i]->lg = d3 * d1 + (1 - d3) * d2;
}
s_x[i]->dg = (c2 * s_x[i]->equiv * s_x[i]->equiv /
((1.0 + s_x[i]->dha * b * muhalf) * (1.0 +
s_x[i]->
dha * b *
muhalf)) +
s_x[i]->dg = coef * (c2 * z * z /
((1.0 + s_x[i]->dha * b * muhalf) * (1.0 + s_x[i]->dha * b * muhalf)) +
s_x[i]->dhb) * LOG_10 * s_x[i]->moles;
}
else if (s_x[i]->exch_gflag == 7 && s_x[i]->alk > 0)
@ -741,17 +742,13 @@ gammas(LDBLE mu)
if (llnl_count_temp > 0)
{
s_x[i]->lg =
-a_llnl * muhalf * s_x[i]->equiv * s_x[i]->equiv /
coef * (-a_llnl * muhalf * z * z /
(1.0 + s_x[i]->dha * b_llnl * muhalf) +
bdot_llnl * mu +
bdot_llnl * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg =
(c2_llnl * s_x[i]->equiv * s_x[i]->equiv /
((1.0 + s_x[i]->dha * b_llnl * muhalf) * (1.0 +
s_x[i]->
dha *
b_llnl *
muhalf)) +
coef * (c2_llnl * z * z /
((1.0 + s_x[i]->dha * b_llnl * muhalf) * (1.0 + s_x[i]->dha * b_llnl * muhalf)) +
bdot_llnl) * LOG_10 * s_x[i]->moles;
}
else
@ -773,16 +770,15 @@ gammas(LDBLE mu)
else
{
if (s_x[i]->alk <= 0)
{
s_x[i]->lg = 0.0;
}
else
{
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
}
s_x[i]->dg = 0.0;
}
}
if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)
gammas_a_f(i); // appt
break;
case 5: /* Always 1.0 */
s_x[i]->lg = 0.0;
@ -879,6 +875,56 @@ gammas(LDBLE mu)
return (OK);
}
/* ------------------------------------------------------------------------------- */
int Phreeqc::gammas_a_f(int i1)
/* ------------------------------------------------------------------------------- */
{
int i, j;
LDBLE d2, d3, coef = 0, sum = 0;
char name[MAX_LENGTH];
//struct master *m_ptr;
i = i1;
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
if (s_x[i]->rxn_x->token[j].s->type == EX)
{
strcpy(name, s_x[i]->rxn_x->token[j].s->name);
//m_ptr = s_x[i]->rxn_x->token[j].s->primary->elt->master; // appt debug
break;
}
}
for (i = 0; i < count_s_x; i++)
{
if (s_x[i]->gflag != 4 || s_x[i]->primary)
continue;
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
if (s_x[i]->rxn_x->token[j].s->type == EX)
{
if (!strcmp(name, s_x[i]->rxn_x->token[j].s->name))
sum += s_x[i]->moles * s_x[i]->equiv;
break;
}
}
}
i = i1;
d2 = s_x[i]->moles * s_x[i]->equiv / sum;
if (d2 > 1) d2 = 1;
//if (iterations > 19)
// i += 0; // appt debug
d3 = 0.5;
if (s_x[i]->a_f > 2)
{
d3 += (s_x[i]->a_f - 2) / 10; if (d3 > 0.8) d3 = 0.8;
}
d2 = s_x[i]->dw_a * d3 + (1 - d3) * d2;
s_x[i]->lg -= s_x[i]->a_f * (1 - d2);
s_x[i]->dw_a = d2;
return 0;
}
/* ------------------------------------------------------------------------------- */
int Phreeqc::
ineq(int in_kode)
/* ------------------------------------------------------------------------------- */
@ -1597,9 +1643,7 @@ ineq(int in_kode)
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
ineq_array[l_count_rows * max_column_count + i] = 1.0;
ineq_array[l_count_rows * max_column_count + count_unknowns] =
//x[i]->moles;
0.99 * x[i]->moles - MIN_TOTAL_SS;
//0.9 * x[i]->moles;
back_eq[l_count_rows] = i;
l_count_rows++;
}
@ -4775,7 +4819,7 @@ initial_guesses(void)
*/
int i;
cxxSolution *solution_ptr;
// mu_x is reset here, but the real, already calculated mu_x must be used for INITIAL_EXCHANGE & _SURFACE appt
solution_ptr = use.Get_solution_ptr();
mu_x =
s_hplus->moles +
@ -5494,11 +5538,11 @@ numerical_jacobian(void)
int i, j;
cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr();
if (!
(numerical_deriv ||
(numerical_deriv ||
(use.Get_surface_ptr() != NULL && use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) ||
(gas_phase_ptr != NULL && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME &&
(gas_phase_ptr->Get_pr_in() || force_numerical_fixed_volume) && numerical_fixed_volume)
))
(gas_phase_ptr != NULL && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME &&
(gas_phase_ptr->Get_pr_in() || force_numerical_fixed_volume) && numerical_fixed_volume)
))
return(OK);
calculating_deriv = TRUE;
@ -5506,9 +5550,9 @@ numerical_jacobian(void)
molalities(TRUE);
mb_sums();
residuals();
/*
* Clear array, note residuals are in array[i, count_unknowns+1]
*/
/*
* Clear array, note residuals are in array[i, count_unknowns+1]
*/
for (i = 0; i < count_unknowns; i++)
{
@ -5517,10 +5561,10 @@ numerical_jacobian(void)
for (i = 1; i < count_unknowns; i++)
{
memcpy((void *) &(my_array[i * (count_unknowns + 1)]),
(void *) &(my_array[0]), (size_t) count_unknowns * sizeof(LDBLE));
(void *) &(my_array[0]), (size_t)count_unknowns * sizeof(LDBLE));
}
base = (LDBLE *) PHRQ_malloc((size_t) count_unknowns * sizeof(LDBLE));
base = (LDBLE *)PHRQ_malloc((size_t)count_unknowns * sizeof(LDBLE));
if (base == NULL)
malloc_error();
for (i = 0; i < count_unknowns; i++)

27
phreeqcpp/pitzer.cpp
View File

@ -1836,7 +1836,7 @@ pitzer_revise_guesses(void)
max_iter = 100;
if (iterations < 0 && (use.Get_surface_in() || use.Get_exchange_in()))
gammas_pz(); // DL_pitz : for SURF estimates
gammas_pz(true); // DL_pitz : for SURF estimates
l_iter = 0;
repeat = TRUE;
fail = FALSE;
@ -1977,7 +1977,6 @@ pitzer_revise_guesses(void)
{
mu_x = 1e-8;
}
//gammas_pz();
return (OK);
}
@ -2077,7 +2076,7 @@ Restart:
d2 = d * mu_x;
mu_x += d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
gammas_pz(false);
break;
case PP:
continue;
@ -2101,7 +2100,7 @@ Restart:
if (max_unknowns > pz_max_unknowns)
{
base = (LDBLE *) free_check_null(base);
gammas_pz();
gammas_pz(false);
jacobian_sums();
goto Restart;
}
@ -2153,7 +2152,7 @@ Restart:
case MU:
mu_x -= d2;
//k_temp(tc_x, patm_x);
gammas(mu_x);
gammas_pz(false);
break;
case GAS_MOLES:
if (gas_in == FALSE)
@ -2285,7 +2284,7 @@ model_pz(void)
/*
* Calculate jacobian
*/
gammas_pz();
gammas_pz(false); // appt: no gammas_a_f here
jacobian_sums();
jacobian_pz();
/*
@ -2314,7 +2313,8 @@ model_pz(void)
}
reset();
}
gammas_pz();
// appt calculate gammas_a_f here
gammas_pz(true);
if (full_pitzer == TRUE)
pitzer();
if (always_full_pitzer == TRUE)
@ -2455,7 +2455,7 @@ check_gammas_pz(void)
/* ---------------------------------------------------------------------- */
int Phreeqc::
gammas_pz()
gammas_pz(bool exch_a_f)
/* ---------------------------------------------------------------------- */
{
/*
@ -2538,7 +2538,7 @@ gammas_pz()
* calculate exchange gammas
*/
if (use.Get_exchange_ptr() != NULL)
if (use.Get_exchange_ptr() != NULL && exch_a_f) // appt for gammas_a_f
{
for (i = 0; i < count_s_x; i++)
{
@ -2560,14 +2560,11 @@ gammas_pz()
* Find CEC
* z contains valence of cation for exchange species, alk contains cec
*/
/* !!!!! */
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
if (s_x[i]->rxn_x->token[j].s->type == EX)
{
s_x[i]->alk =
s_x[i]->rxn_x->token[j].s->primary->unknown->
moles;
s_x[i]->alk = s_x[i]->rxn_x->token[j].s->primary->unknown->moles;
break;
}
}
@ -2598,9 +2595,11 @@ gammas_pz()
continue;
coef = s_x[i]->rxn_x->token[j].coef;
s_x[i]->lg += coef * s_x[i]->rxn_x->token[j].s->lg;
s_x[i]->dg += coef * s_x[i]->rxn_x->token[j].s->dg;
/*s_x[i]->dg += coef * s_x[i]->rxn_x->token[j].s->dg;*//* remove? */
}
}
if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)
gammas_a_f(i); // appt
}
}
}

10
phreeqcpp/prep.cpp
View File

@ -1405,11 +1405,11 @@ build_model(void)
write_mass_action_eqn_x(STOP);
trxn_reverse_k();
rxn_free(phases[i]->rxn_x);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\nPhase: %s\n", phases[i]->name));
trxn_print();
}
//if (debug_prep == TRUE)
//{
// output_msg(sformatf( "\nPhase: %s\n", phases[i]->name));
// trxn_print();
//}
phases[i]->rxn_x = rxn_alloc(count_trxn + 1);
trxn_copy(phases[i]->rxn_x);
write_phase_sys_total(i);

3
phreeqcpp/print.cpp
View File

@ -520,8 +520,7 @@ print_exchange(void)
/*
* Print species data
*/
/* !!!!! */
if (master_ptr->total > 1.0e-10)
if (master_ptr->total > 1.0e-16)
{
if (species_list[i].s->equiv != 0.0)
{

4
phreeqcpp/read.cpp
View File

@ -45,7 +45,7 @@ read_input(void)
Rxn_new_ss_assemblage.clear();
Rxn_new_surface.clear();
Rxn_new_temperature.clear(); // not used
phrq_io->Set_echo_on(true); // **appt
phrq_io->Set_echo_on(true);
/*
* Initialize keyword counters
*/
@ -8715,7 +8715,7 @@ read_print(void)
if (pr.echo_input == 0)
phrq_io->Set_echo_on(false);
else
phrq_io->Set_echo_on(true); //**appt
phrq_io->Set_echo_on(true);
break;
case 32: /* warning */
case 33: /* warnings */

10
phreeqcpp/step.cpp
View File

@ -1102,6 +1102,7 @@ add_ss_assemblage(cxxSSassemblage *ss_assemblage_ptr)
{
char * token = string_duplicate(phase_ptr->formula);
ptr = &(token[0]);
count_elts = 0; // appt
get_elts_in_species(&ptr, 1.0);
free_check_null(token);
for (k = 0; k < count_elts; k++)
@ -1470,10 +1471,13 @@ solution_check(void)
"Element %s has negative moles in solution, %e. \n\tErroneous mole balance occurs as moles are added to produce zero moles.\n\tUsually caused by KINETICS, REACTION, or diffuse layer calculation.\n\tMay be due to large time steps in early part of KINETICS simulation or negative concentrations in the diffuse layer.",
master_ptr->elt->name, (LDBLE) master_ptr->total);
*/
error_string = sformatf(
if (state != TRANSPORT)
{
error_string = sformatf(
"Negative moles in solution %d for %s, %e. Recovering...",
cell_no, master_ptr->elt->name, (double) master_ptr->total);
warning_msg(error_string);
cell_no, master_ptr->elt->name, (double)master_ptr->total);
warning_msg(error_string);
}
return (MASS_BALANCE);
}
return (OK);

10
phreeqcpp/transport.cpp
View File

@ -1,4 +1,4 @@
#include "Utils.h"
#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "Exchange.h"
@ -1111,7 +1111,7 @@ init_mix(void)
if (2.25 * maxmix + 1.0 > (double)INT_MAX)
{
char token[MAX_LENGTH];
sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step.", 2.25 * maxmix);
sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step, or increase cell-lengths.", 2.25 * maxmix);
error_msg(token, STOP);
}
if (bcon_first == 1 || bcon_last == 1)
@ -1227,6 +1227,12 @@ init_mix(void)
l_nmix = 0;
else
{
if (1.5 * maxmix > (double)INT_MAX)
{
char token[MAX_LENGTH];
sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step, or increase cell-lengths.", 1.5 * maxmix);
error_msg(token, STOP);
}
l_nmix = 1 + (int) floor(1.5 * maxmix);
if ((ishift != 0) && ((bcon_first == 1) || (bcon_last == 1)))