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

Merge commit 'd4dfe590eb539c35219174adf5314480eebdf7f8'

Browse Source
This commit is contained in:
Darth Vader 2021-08-17 15:36:31 +00:00
commit 5752c54f20
14 changed files with 627 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/phreeqcpp/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
src/phreeqcpp/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);

28
src/phreeqcpp/gases.cpp
View File

@ -101,7 +101,7 @@ build_fixed_volume_gas(void)
}
/* All elements in gas */
for (j = 0; j < count_elts; j++)
for (j = 0; j < (int) count_elts; j++)
{
unknown_ptr = NULL;
if (strcmp(elt_list[j].elt->name, "H") == 0)
@ -149,7 +149,7 @@ build_fixed_volume_gas(void)
output_msg(sformatf( "\n\tJacobian summations %s.\n\n",
phase_ptr->name));
}
for (j = 0; j < count_elts; j++)
for (j = 0; j < (int) count_elts; j++)
{
unknown_ptr = NULL;
if (strcmp(elt_list[j].elt->name, "H") == 0)
@ -436,23 +436,31 @@ 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)") || !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;
else if (!strcmp(phase_ptr1->name, "Propane(g)"))
a_aa *= 0.45;
}
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)") || !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;
else if (!strcmp(phase_ptr->name, "Propane(g)"))
a_aa *= 0.45;
}
a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa;
a_aa_sum2 += phase_ptr1->fraction_x * a_aa;
@ -598,10 +606,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
src/phreeqcpp/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
src/phreeqcpp/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
src/phreeqcpp/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)

274
src/phreeqcpp/model.cpp
View File

@ -2613,13 +2613,12 @@ calc_gas_pressures(void)
V_m = (2. * gas_phase_ptr->Get_v_m() + V_m) / 3;
else
V_m = (1. * gas_phase_ptr->Get_v_m() + V_m) / 2;
if (iterations > 99 && numerical_fixed_volume == false)
if ((pitzer_model || iterations > 99) && numerical_fixed_volume == false)
{
//V_m *= 1; /* debug */
numerical_fixed_volume = true;
//switch_numerical = true;
warning_msg
("Numerical method failed, switching to numerical derivatives.");
if (!pitzer_model) warning_msg("Numerical method failed, switching to numerical derivatives.");
prep();
//switch_numerical = false;
}
@ -2649,8 +2648,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)
{
@ -5305,7 +5304,8 @@ ss_f(LDBLE xb, LDBLE l_a0, LDBLE l_a1, LDBLE l_kc, LDBLE l_kb, LDBLE xcaq,
f = xcaq * (xb / r + xc) + xbaq * (xb + r * xc) - 1;
return (f);
}
//#define ORIGINAL
#ifdef ORIGINAL
/* ---------------------------------------------------------------------- */
int Phreeqc::
numerical_jacobian(void)
@ -5425,8 +5425,8 @@ numerical_jacobian(void)
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
d2 = (x[i]->moles > 1 ? 1 : 20);
d2 *= d * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
x[i]->moles += d2;
@ -5546,7 +5546,265 @@ numerical_jacobian(void)
calculating_deriv = FALSE;
return OK;
}
#else
/* ---------------------------------------------------------------------- */
int Phreeqc::
numerical_jacobian(void)
/* ---------------------------------------------------------------------- */
{
std::vector<double> base;
LDBLE d, d1, d2;
int i, j;
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
std::vector<class phase*> phase_ptrs;
std::vector<class phase> base_phases;
double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
if (!
(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)
))
return(OK);
//jacobian_sums();
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
//cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
base_gas_phase = *gas_phase_ptr;
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
{
const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
phase_ptrs.push_back(phase_ptr);
base_phases[i] = *phase_ptr;
}
}
calculating_deriv = TRUE;
gammas(mu_x);
molalities(TRUE);
mb_sums();
//mb_gases();
//mb_ss();
residuals();
/*
* Clear array, note residuals are in array[i, count_unknowns+1]
*/
//for (i = 0; i < count_unknowns; i++)
//{
// my_array[i] = 0.0;
//}
//for (i = 1; i < count_unknowns; i++)
//{
// memcpy((void*)&(my_array[(size_t)i * (count_unknowns + 1)]),
// (void*)&(my_array[0]), count_unknowns * sizeof(LDBLE));
//}
base.resize(count_unknowns);
base = residual;
d = 0.0001;
d1 = d * LOG_10;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
{
switch (x[i]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
x[i]->master[0]->s->la += d;
d2 = d;// *LOG_10;
break;
case MH:
s_eminus->la += d;
d2 = d;// *LOG_10;
break;
case AH2O:
x[i]->master[0]->s->la += d;
d2 = d;// *LOG_10;
break;
case PITZER_GAMMA:
x[i]->s->lg += d;
d2 = d;
break;
case MH2O:
//mass_water_aq_x *= (1 + d);
//x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
//d2 = log(1.0 + d);
//break;
// DL_pitz
d1 = mass_water_aq_x * d;
mass_water_aq_x += d1;
if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL)
mass_water_bulk_x += d1;
x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
d2 = d1;
break;
case MU:
d2 = d * mu_x;
mu_x += d2;
//gammas(mu_x);
break;
case PP:
for (j = 0; j < count_unknowns; j++)
{
delta[j] = 0.0;
}
d2 = -1e-8;
delta[i] = d2;
reset();
d2 = delta[i];
break;
case SS_MOLES:
if (x[i]->ss_in == FALSE)
continue;
for (j = 0; j < count_unknowns; j++)
{
delta[j] = 0.0;
}
/*d2 = -1e-8; */
d2 = d * 10 * x[i]->moles;
//d2 = -.1 * x[i]->moles;
/*
if (d2 > -1e-10) d2 = -1e-10;
calculating_deriv = FALSE;
*/
delta[i] = d2;
/*fprintf (stderr, "delta before reset %e\n", delta[i]); */
reset();
d2 = delta[i];
/*fprintf (stderr, "delta after reset %e\n", delta[i]); */
break;
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = (x[i]->moles > 1 ? 1 : 30);
d2 *= d * x[i]->moles;
d2 = (d2 < ineq_tol ? ineq_tol : d2);
//if (d2 < 1e-14)
// d2 = 1e-14;
x[i]->moles += d2;
break;
}
gammas(mu_x);
molalities(TRUE);
mb_sums();
//mb_gases();
//mb_ss();
residuals();
for (j = 0; j < count_unknowns; j++)
{
my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] = -(residual[j] - base[j]) / d2;
if (x[i]->type == MH2O) // DL_pitz
my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] *= mass_water_aq_x;
}
switch (x[i]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
case AH2O:
x[i]->master[0]->s->la -= d2;
break;
case MH:
s_eminus->la -= d2;
if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0)
{
/*output_msg(sformatf( "Zero diagonal for MH\n")); */
my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] =
under(s_h2->lm) * 2;
}
break;
case PITZER_GAMMA:
x[i]->s->lg -= d2;
break;
case MH2O:
//mass_water_aq_x /= (1 + d);
//x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
//break;
//DL_pitz
mass_water_aq_x -= d2;
if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL)
mass_water_bulk_x -= d2;
x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
break;
case MU:
mu_x -= d2;
//gammas(mu_x);
break;
case PP:
delta[i] = -d2;
reset();
break;
case SS_MOLES:
delta[i] = -d2;
reset();
break;
case GAS_MOLES:
x[i]->moles -= d2;
break;
}
if (use.Get_surface_ptr() != NULL)
{
*use.Get_surface_ptr() = base_surface;
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;
for (size_t g = 0; g < base_phases.size(); g++)
{
*phase_ptrs[g] = base_phases[g];
}
}
//gammas(mu_x);
//molalities(TRUE);
//mb_sums();
////mb_gases();
////mb_ss();
//residuals();
}
gammas(mu_x);
molalities(TRUE);
mb_sums();
//mb_gases();
//mb_ss();
residuals();
//for (i = 0; i < count_unknowns; i++)
//{
// //Debugging
// if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
// fabs(residual[i]) + fabs(base[i]) > 1e-8)
// {
// std::cerr << iterations << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl;
// }
//}
base.clear();
calculating_deriv = FALSE;
return OK;
}
#endif
/* ---------------------------------------------------------------------- */
void Phreeqc::
set_inert_moles(void)

67
src/phreeqcpp/pitzer.cpp
View File

@ -1948,19 +1948,40 @@ jacobian_pz(void)
/* ---------------------------------------------------------------------- */
{ // calculate the derivatives numerically
std::vector<double> base;
std::vector<class phase*> phase_ptrs;
std::vector<class phase> base_phases;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
LDBLE d, d1, d2;
int i, j;
calculating_deriv = 1;
Restart:
size_t pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x);
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
base_gas_phase = *gas_phase_ptr;
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
{
const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
phase_ptrs.push_back(phase_ptr);
base_phases[i] = *phase_ptr;
}
}
calculating_deriv = 1;
molalities(TRUE);
if (full_pitzer == TRUE)
{
molalities(TRUE);
pitzer();
residuals();
}
mb_sums();
residuals();
size_t pz_max_unknowns = max_unknowns;
base.resize(count_unknowns);
for (i = 0; i < count_unknowns; i++)
{
@ -2026,9 +2047,11 @@ Restart:
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
d2 = (x[i]->moles > 1 ? 1 : 30);
d2 *= d * x[i]->moles;
d2 = (d2 < ineq_tol ? ineq_tol : d2);
//if (d2 < 1e-14)
// d2 = 1e-14;
x[i]->moles += d2;
break;
case MU:
@ -2124,12 +2147,38 @@ Restart:
reset();
break;
}
if (use.Get_surface_ptr() != NULL)
{
*use.Get_surface_ptr() = base_surface;
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;
for (size_t g = 0; g < base_phases.size(); g++)
{
*phase_ptrs[g] = base_phases[g];
}
}
//molalities(TRUE);
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
}
molalities(TRUE);
if (full_pitzer == TRUE)
pitzer();
mb_sums();
residuals();
//for (i = 0; i < count_unknowns; i++)
//{
// //Debugging
// if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
// fabs(residual[i]) + fabs(base[i]) > 1e-6)
// {
// std::cerr << iterations << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl;
// }
//}
base.clear();
calculating_deriv = 0;
return OK;

34
src/phreeqcpp/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
src/phreeqcpp/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
src/phreeqcpp/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;
}

210
src/phreeqcpp/sit.cpp
View File

@ -798,7 +798,8 @@ sit_revise_guesses(void)
/*gammas(mu_x); */
return (OK);
}
//#define ORIGINAL
#ifdef ORIGINAL
/* ---------------------------------------------------------------------- */
int Phreeqc::
jacobian_sit(void)
@ -870,7 +871,8 @@ Restart:
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = d * x[i]->moles;
d2 = (x[i]->moles > 1 ? 1 : 20);
d2 *= d * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
x[i]->moles += d2;
@ -951,7 +953,209 @@ Restart:
residuals();
return OK;
}
#else
/* ---------------------------------------------------------------------- */
int Phreeqc::
jacobian_sit(void)
/* ---------------------------------------------------------------------- */
{
std::vector<double> base;
LDBLE d, d1, d2;
int i, j;
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
std::vector<class phase*> phase_ptrs;
std::vector<class phase> base_phases;
double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
cxxGasPhase base_gas_phase;
cxxSurface base_surface;
Restart:
if (use.Get_surface_ptr() != NULL)
{
base_surface = *use.Get_surface_ptr();
}
if (use.Get_gas_phase_ptr() != NULL)
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
base_gas_phase = *gas_phase_ptr;
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
{
const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
phase_ptrs.push_back(phase_ptr);
base_phases[i] = *phase_ptr;
}
}
calculating_deriv = 1;
size_t pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x);
molalities(TRUE);
if (full_pitzer == TRUE)
{
sit();
}
mb_sums();
residuals();
base = residual; // std::vectors
d = 0.0001;
d1 = d * LOG_10;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
{
switch (x[i]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
x[i]->master[0]->s->la += d;
d2 = d1;
break;
case AH2O:
x[i]->master[0]->s->la += d;
d2 = d1;
break;
case PITZER_GAMMA:
if (!full_pitzer)
continue;
x[i]->s->lg += d;
d2 = d;
break;
case MH2O:
mass_water_aq_x *= (1.0 + d);
x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
d2 = log(1.0 + d);
break;
case MH:
s_eminus->la += d;
d2 = d1;
break;
/*
if (pitzer_pe == TRUE)
{
s_eminus->la += d;
d2 = d1;
break;
}
else
{
continue;
}
*/
case GAS_MOLES:
if (gas_in == FALSE)
continue;
d2 = (x[i]->moles > 1 ? 1 : 20);
d2 *= d * x[i]->moles;
if (d2 < 1e-14)
d2 = 1e-14;
x[i]->moles += d2;
break;
case MU:
//continue;
d2 = d * mu_x;
mu_x += d2;
//k_temp(tc_x, patm_x);
gammas_sit();
break;
case PP:
case SS_MOLES:
continue;
break;
}
molalities(TRUE);
if (max_unknowns > pz_max_unknowns)
{
gammas_sit();
jacobian_sums();
goto Restart;
}
if (full_pitzer == TRUE)
sit();
mb_sums();
residuals();
for (j = 0; j < count_unknowns; j++)
{
my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] =
-(residual[j] - base[j]) / d2;
}
switch (x[i]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
case AH2O:
x[i]->master[0]->s->la -= d;
break;
case MH:
s_eminus->la -= d;
if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0)
{
my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] =
exp(s_h2->lm * LOG_10) * 2;
}
break;
case PITZER_GAMMA:
x[i]->s->lg -= d;
break;
case MH2O:
mass_water_aq_x /= (1 + d);
x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
break;
case MU:
mu_x -= d2;
//k_temp(tc_x, patm_x);
gammas_sit();
break;
case GAS_MOLES:
if (gas_in == FALSE)
continue;
x[i]->moles -= d2;
break;
}
if (use.Get_surface_ptr() != NULL)
{
*use.Get_surface_ptr() = base_surface;
}
if (use.Get_gas_phase_ptr() != NULL)
{
*use.Get_gas_phase_ptr() = base_gas_phase;
for (size_t g = 0; g < base_phases.size(); g++)
{
*phase_ptrs[g] = base_phases[g];
}
}
}
molalities(TRUE);
if (full_pitzer == TRUE)
sit();
mb_sums();
residuals();
//for (i = 0; i < count_unknowns; i++)
//{
// //Debugging
// if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
// fabs(residual[i]) + fabs(base[i]) > 1e-8)
// {
// std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl;
// }
//}
calculating_deriv = 0;
return OK;
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
model_sit(void)

12
src/phreeqcpp/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
src/phreeqcpp/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;