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

Tony H2S. Amm.dat, phreeqc.dat, pitzer.dat, utf8, updated test cases

Browse Source
This commit is contained in:
David Parkhurst 2021-07-06 21:52:06 -06:00
parent 48cb5e8969
commit 9022ded877
14 changed files with 94 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
Solution.cxx
View File

@ -114,6 +114,17 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
this->n_user = this->n_user_end = l_n_user;
this->new_def = false;
this->ah2o = 0;
// potV is an external variable, imposed in a given solution, not mixed.
std::map < int, cxxSolution >::const_iterator sol = solutions.find(mix.Get_n_user());
const cxxSolution *cxxsoln_ptr1;
if (sol != solutions.end())
{
cxxsoln_ptr1 = &(sol->second);
if (cxxsoln_ptr1->new_def)
this->potV = 0.0;
else
this->potV = cxxsoln_ptr1->potV;
}
//
// Mix solutions
//
@ -121,8 +132,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
std::map < int, LDBLE >::const_iterator it;
for (it = mixcomps.begin(); it != mixcomps.end(); it++)
{
std::map < int, cxxSolution >::const_iterator sol =
solutions.find(it->first);
sol = solutions.find(it->first);
if (sol == solutions.end())
{
std::ostringstream msg;
@ -131,7 +141,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
}
else
{
const cxxSolution *cxxsoln_ptr1 = &(sol->second);
cxxsoln_ptr1 = &(sol->second);
this->add(*cxxsoln_ptr1, it->second);
}
}
@ -1388,7 +1398,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
this->cb += addee.cb * extensive;
this->density = f1 * this->density + f2 * addee.density;
this->patm = f1 * this->patm + f2 * addee.patm;
this->potV = f1 * this->potV + f2 * addee.potV;
// 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;

14
class_main.cpp
View File

@ -147,8 +147,8 @@ main_method(int argc, char *argv[])
}
Phreeqc MyCopy = *this;
this->clean_up();
//this->init();
//this->initialize();
// this->init();
// this->initialize();
/*
* Read input data for simulation
*/
@ -167,7 +167,7 @@ main_method(int argc, char *argv[])
* Display successful status
*/
pr.headings = TRUE;
//errors = do_status();
// errors = do_status();
if (errors != 0)
{
return errors;
@ -292,7 +292,7 @@ write_banner(void)
/* version */
#ifdef NPP
len = sprintf(buffer, "* PHREEQC-%s *", "3.6.5");
len = sprintf(buffer, "* PHREEQC-%s *", "3.7.1");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -316,7 +316,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "February 24, 2021");
len = sprintf(buffer, "%s", "July 5, 2021");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -329,6 +329,7 @@ write_banner(void)
return 0;
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
process_file_names(int argc, char *argv[], std::istream **db_cookie,
@ -499,7 +500,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
output_msg(sformatf(" Input file: %s\n", in_file.c_str()));
output_msg(sformatf(" Output file: %s\n", out_file.c_str()));
#ifdef NPP
output_msg(sformatf("Using PHREEQC: version 3.6.5, compiled February 24, 2021\n"));
output_msg(sformatf("Using PHREEQC: version 3.7.1, compiled July 5, 2021\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token.c_str()));
#ifdef NPP
@ -508,6 +509,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
/*
* local cleanup
*/
line = (char *) free_check_null(line);
line_save = (char *) free_check_null(line_save);

16
gases.cpp
View File

@ -436,22 +436,22 @@ calc_PR(void)
{
if (!strcmp(phase_ptr1->name, "CO2(g)"))
a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217
else if (!strcmp(phase_ptr1->name, "H2S(g)"))
else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr1->name, "CH4(g)"))
else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "N2(g)"))
else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)"))
a_aa *= 0.51;
}
if (!strcmp(phase_ptr1->name, "H2O(g)"))
{
if (!strcmp(phase_ptr->name, "CO2(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "H2S(g)"))
else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "CH4(g)"))
else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "N2(g)"))
else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)"))
a_aa *= 0.51;
}
a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa;
@ -598,10 +598,10 @@ calc_PR(void)
{
phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
}
else
phi = -3.0; // fugacity coefficient > 0.05
phi = -3.0; // fugacity coefficient = 0.05
phase_ptr->pr_phi = exp(phi);
phase_ptr->pr_si_f = phi / LOG_10; // pr_si_f updated
// ****

1
global_structures.h
View File

@ -41,6 +41,7 @@
#define DISP 2
#define STAG 3
#define NOMIX 4
#define MIX_BS 5 // mix boundary solutions in electromigration
#define CONVERGED 2
#define MASS_BALANCE 3

4
inverse.cpp
View File

@ -272,7 +272,7 @@ setup_inverse(class inverse *inv_ptr)
/*
* Define inv_zero and inv_zero array, delta
*/
for (i = 0; i < max; i++)
for (i = 0; i < (int) max; i++)
inv_zero[i] = 0.0;
memcpy((void *) &(delta[0]), (void *) &(inv_zero[0]),
@ -3595,7 +3595,7 @@ check_isotopes(class inverse *inv_ptr)
/* use solution-defined uncertainties second */
}
#ifdef NPP
else if (inv_ptr->i_u[i].count_uncertainties > 0
else if (inv_ptr->i_u[i].uncertainties.size() > 0
&& !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
#else
else if (inv_ptr->i_u[i].uncertainties.size() > 0

4
kinetics.cpp
View File

@ -1595,7 +1595,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver)
/*
* i --user number for soln, reaction, etc.
* use_mix --integer flag
state == TRANSPORT: DISP, STAG, NOMIX
state == TRANSPORT: DISP, STAG, NOMIX, MIX_BS
state == REACTION: TRUE, FALSE
* use_kinetics --true or false flag to calculate kinetic reactions
* nsaver --user number to store solution
@ -1615,7 +1615,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver)
use.Set_n_mix_user(i);
use.Set_n_mix_user_orig(i);
}
else if (use_mix == STAG && multi_Dflag != TRUE)
else if ((use_mix == STAG && multi_Dflag != TRUE) || use_mix == MIX_BS)
{
use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i));
if (use.Get_mix_ptr() != NULL)

7
model.cpp
View File

@ -2649,8 +2649,8 @@ calc_gas_pressures(void)
lp += rxn_ptr->s->la * rxn_ptr->coef;
}
phase_ptr->p_soln_x = exp(LOG_10 * (lp - phase_ptr->pr_si_f));
if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90)
phase_ptr->p_soln_x = 90;
//if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90)
// phase_ptr->p_soln_x = 90;
if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE)
{
@ -5425,8 +5425,7 @@ numerical_jacobian(void)
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
d2 = d * 20 * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
x[i]->moles += d2;

7
pitzer.cpp
View File

@ -2026,9 +2026,10 @@ Restart:
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
d2 = d * 30 * x[i]->moles;
d2 = (d2 < ineq_tol ? ineq_tol : d2);
//if (d2 < 1e-14)
// d2 = 1e-14;
x[i]->moles += d2;
break;
case MU:

34
prep.cpp
View File

@ -3841,15 +3841,11 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
{
if (!strcmp(phase_ptr1->name, "CO2(g)"))
a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217
else if (!strcmp(phase_ptr1->name, "H2S(g)"))
else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr1->name, "CH4(g)"))
else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "Mtg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "Methane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "N2(g)"))
else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "Ethane(g)"))
a_aa *= 0.51;
@ -3860,15 +3856,11 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
{
if (!strcmp(phase_ptr->name, "CO2(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "H2S(g)"))
else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)"))
a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "CH4(g)"))
else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "Mtg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "Methane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "N2(g)"))
else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "Ethane(g)"))
a_aa *= 0.51;
@ -4002,17 +3994,17 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
{
phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
//if (phi > 4.44)
// phi = 4.44;
}
else
phi = -3.0; // fugacity coefficient > 0.05
if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3)
{
// avoid such phi...
phi = -3;
}
phi = -3.0; // fugacity coefficient = 0.05
//if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3)
//{
//// avoid such phi...
// phi = -3;
//}
phase_ptr->pr_phi = exp(phi);
phase_ptr->pr_si_f = phi / LOG_10;
// for initial equilibrations, adapt log_k of the gas phase...

4
print.cpp
View File

@ -708,8 +708,8 @@ print_gas_phase(void)
(double) initial_moles,
(double) moles,
(double) delta_moles));
if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90)
output_msg(" WARNING: The pressure of H2O(g) is above the program limit: use the polynomial for log_k.\n");
//if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90)
// output_msg(" WARNING: The pressure of H2O(g) is fixed to the program limit.\n");
}
output_msg("\n");

14
read.cpp
View File

@ -2619,7 +2619,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
/*
* Read supcrt parms and Ionic strength terms
*/
for (j = 0; j < 9; j++)
for (j = 0; j < 11; j++)
{
delta_v[j] = 0.0;
}
@ -2627,7 +2627,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
/* Vmax, dmax...
delta_v[10] = 999.0;
delta_v[11] = 1.0; */
j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT /*SCANFORMAT SCANFORMAT */,
j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT ,
/* a1..a4 */
&(delta_v[0]), &(delta_v[1]), &(delta_v[2]), &(delta_v[3]),
/* wref */
@ -2635,9 +2635,9 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
/* b_Av */
&(delta_v[5]),
/* c1..c4 */
&(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]));
/* vmax, dmax
&(delta_v[10]), &(delta_v[11])); */
&(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]),
//vmP, exP
&(delta_v[10]), &(delta_v[11]));
if (j < 1)
{
input_error++;
@ -2665,8 +2665,8 @@ read_vm_only(const char* cptr, LDBLE * delta_v, DELTA_V_UNIT * units)
char token[MAX_LENGTH];
/*
* Read analytical expression
*/
for (j = 0; j < 8; j++)
*/
for (j = 0; j < 9; j++)
{
delta_v[j] = 0.0;
}

2
sit.cpp
View File

@ -870,7 +870,7 @@ Restart:
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
d2 = d * 20 * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
x[i]->moles += d2;

12
step.cpp
View File

@ -57,11 +57,19 @@ step(LDBLE step_fraction)
*/
if (use.Get_mix_ptr() != NULL)
{
add_mix(use.Get_mix_ptr());
add_mix(use.Get_mix_ptr());
int n = use.Get_n_mix_user_orig();
if (n == 0 || n == count_cells + 1)
{
cxxSolution *solution_ptr = Utilities::Rxn_find(Rxn_solution_map, n);
if (solution_ptr != NULL && !solution_ptr->Get_new_def())
potV_x = solution_ptr->Get_potV();
}
}
else if (use.Get_solution_ptr() != NULL)
{
add_solution(use.Get_solution_ptr(), 1.0, 1.0);
potV_x = use.Get_solution_ptr()->Get_potV();
cell_no = use.Get_n_solution_user();
}
else
@ -367,7 +375,7 @@ add_solution(cxxSolution *solution_ptr, LDBLE extensive, LDBLE intensive)
tc_x += solution_ptr->Get_tc() * intensive;
ph_x += solution_ptr->Get_ph() * intensive;
patm_x += solution_ptr->Get_patm() * intensive;
potV_x += solution_ptr->Get_potV() * intensive;
//potV_x += solution_ptr->Get_potV() * intensive;
solution_pe_x += solution_ptr->Get_pe() * intensive;
mu_x += solution_ptr->Get_mu() * intensive;
ah2o_x += solution_ptr->Get_ah2o() * intensive;

23
transport.cpp
View File

@ -1,4 +1,4 @@
#include "Utils.h"
#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "Exchange.h"
@ -605,7 +605,10 @@ transport(void)
if (i == 0 || i == count_cells + 1)
{
run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
if (dV_dcell)
run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
else
run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
}
else
{
@ -868,7 +871,12 @@ transport(void)
status(0, token);
if (i == 0 || i == count_cells + 1)
run_reactions(i, kin_time, NOMIX, step_fraction);
{
if (dV_dcell)
run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
else
run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
}
else
{
run_reactions(i, kin_time, DISP, step_fraction);
@ -1070,6 +1078,12 @@ print_punch(int i, boolean active)
if (!active)
run_reactions(i, 0, NOMIX, 0);
cell_no = i;
if (dV_dcell)
{
use.Set_n_solution_user(i);
use.Get_solution_ptr()->Set_potV(cell_data[i].potV);
potV_x = cell_data[i].potV;
}
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, i));
if (use.Get_kinetics_ptr() != NULL)
{
@ -1413,7 +1427,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction)
for (std::map < int, LDBLE >::const_iterator it = use.Get_mix_ptr()->Get_mixComps().begin();
it != use.Get_mix_ptr()->Get_mixComps().end(); it++)
{
if (it->first > i && it->first < all_cells)
if (it->first > i && it->first < all_cells && it->first != count_cells + 1)
{
k = it->first;
ptr_imm = Utilities::Rxn_find(Rxn_solution_map, k);
@ -6053,6 +6067,7 @@ calc_vm_Cl(void)
V_Cl = s_ptr->logk[vma1] + s_ptr->logk[vma2] / pb_s +
(s_ptr->logk[vma3] + s_ptr->logk[vma4] / pb_s) / TK_s -
s_ptr->logk[wref] * QBrn;
/* the ionic strength term * I^0.5... */
if (s_ptr->logk[b_Av] < 1e-5)
V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu;