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

Squashed 'phreeqcpp/' changes from 97a0327..e7c12e1

Browse Source 
e7c12e1 Merge commit '823f839bc515d692be4dba290f43673b3c0493ab'
823f839 Squashed 'common/' changes from 04cb71a..6e625e5
c0c3f59 Merge remote-tracking branch 'github/master'
6a30753 changed math.h to cmath.
171ff7f ifdef for NAN
f0839d1 fixed warnings from Intel C++ 2023 compiler
49d0a56 fixed surface.cxx correct_gc. tweaked P_Vm, but it is not too stable.
f2b8caa added viscosity to solution.h, solution.cpp, and xsolution_save
7ff51aa Tony's viscosity with many examples

git-subtree-dir: phreeqcpp
git-subtree-split: e7c12e191414a87ddcedcf73f17d7765702cb75b
This commit is contained in:
Darth Vader 2023-04-10 18:27:31 +00:00
parent 0a7874068e
commit aff166b184
30 changed files with 614 additions and 306 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ChartObject.cpp
View File

@ -11,7 +11,8 @@
#include "ChartObject.h"
#include "Parser.h"
#include <fstream>
#include <math.h>
//#include <math.h>
#include <cmath>
#include <iomanip>
#include "phqalloc.h"

2
PBasic.cpp
View File

@ -7215,7 +7215,7 @@ _NilCheck(void)
return _Escape(-3);
}
#ifdef SKIP
#ifdef NPP
/* The following is suitable for the HP Pascal operating system.
It might want to be revised when emulating another system. */

3
PBasic.h
View File

@ -7,7 +7,8 @@
#include <stdio.h>
#include <limits.h>
#include <ctype.h>
#include <math.h>
//#include <math.h>
#include <cmath>
#include <setjmp.h>
#include "phrqtype.h"
#include "PHRQ_base.h"

4
Phreeqc.h
View File

@ -283,12 +283,12 @@ public:
int sum_diffuse_layer(cxxSurfaceCharge* surface_charge_ptr1);
int calc_all_donnan(void);
int calc_init_donnan(void);
LDBLE calc_psi_avg(cxxSurfaceCharge * charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr);
LDBLE g_function(LDBLE x_value);
LDBLE midpnt(LDBLE x1, LDBLE x2, int n);
void polint(LDBLE* xa, LDBLE* ya, int n, LDBLE xv, LDBLE* yv,
LDBLE* dy);
LDBLE qromb_midpnt(cxxSurfaceCharge* charge_ptr, LDBLE x1, LDBLE x2);
LDBLE calc_psi_avg(cxxSurfaceCharge* charge_ptr, LDBLE surf_chrg_eq);
// inverse.cpp -------------------------------
int inverse_models(void);
@ -1851,7 +1851,7 @@ isfinite handling
#elif defined(HAVE_FINITE)
# define PHR_ISFINITE(x) finite(x)
#elif defined(HAVE_ISNAN)
# define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!isnan(x)) && ((x) != (2.0 * (x)))) )
# define PHR_ISFINITE(x) ( ((x) == 0.0) || ((!std::isnan(x)) && ((x) != (2.0 * (x)))) )
#else
# define PHR_ISFINITE(x) ( ((x) == 0.0) || (((x) == (x)) && ((x) != (2.0 * (x)))) )
#endif

23
Solution.cxx
View File

@ -47,6 +47,7 @@ cxxSolution::cxxSolution(PHRQ_io * io)
this->total_o = 55.55;
this->cb = 0.0;
this->density = 1.0;
this->viscosity = 1.0;
this->mass_water = 1.0;
this->soln_vol = 1.0;
this->total_alkalinity = 0.0;
@ -80,6 +81,7 @@ cxxSolution::operator =(const cxxSolution &rhs)
this->total_h = rhs.total_h;
this->total_o = rhs.total_o;
this->density = rhs.density;
this->viscosity = rhs.viscosity;
this->cb = rhs.cb;
this->mass_water = rhs.mass_water;
this->soln_vol = rhs.soln_vol;
@ -269,6 +271,10 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
s_oss << indent1;
s_oss << "-density " << this->density << "\n";
// new identifier
s_oss << indent1;
s_oss << "-viscosity " << this->viscosity << "\n";
// soln_total conc structures
s_oss << indent1;
s_oss << "-totals" << "\n";
@ -1070,6 +1076,16 @@ cxxSolution::read_raw(CParser & parser, bool check)
opt_save = 27;
}
break;
case 28: // viscosity
if (!(parser.get_iss() >> this->viscosity))
{
this->viscosity = 1.0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for viscosity.",
PHRQ_io::OT_CONTINUE);
}
opt_save = CParser::OPT_DEFAULT;
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
@ -1365,6 +1381,7 @@ cxxSolution::zero()
this->total_o = 0.0;
this->cb = 0.0;
this->density = 1.0;
this->viscosity = 1.0;
this->mass_water = 0.0;
this->soln_vol = 0.0;
this->total_alkalinity = 0.0;
@ -1397,6 +1414,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
this->total_o += addee.total_o * extensive;
this->cb += addee.cb * extensive;
this->density = f1 * this->density + f2 * addee.density;
this->viscosity = f1 * this->viscosity + f2 * addee.viscosity;
this->patm = f1 * this->patm + f2 * addee.patm;
// this->potV = f1 * this->potV + f2 * addee.potV; // appt
this->mass_water += addee.mass_water * extensive;
@ -1574,6 +1592,7 @@ cxxSolution::Serialize(Dictionary & dictionary, std::vector < int >&ints,
doubles.push_back(this->cb);
doubles.push_back(this->mass_water);
doubles.push_back(this->density);
doubles.push_back(this->viscosity);
doubles.push_back(this->soln_vol);
doubles.push_back(this->total_alkalinity);
/*
@ -1660,6 +1679,7 @@ cxxSolution::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
this->cb = doubles[dd++];
this->mass_water = doubles[dd++];
this->density = doubles[dd++];
this->viscosity = doubles[dd++];
this->soln_vol = doubles[dd++];
this->total_alkalinity = doubles[dd++];
/*
@ -1752,6 +1772,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
std::vector< std::string >::value_type("species_map"), // 24
std::vector< std::string >::value_type("log_gamma_map"), // 25
std::vector< std::string >::value_type("potential"), // 26
std::vector< std::string >::value_type("log_molalities_map") // 27
std::vector< std::string >::value_type("log_molalities_map"), // 27
std::vector< std::string >::value_type("viscosity") // 28
};
const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

3
Solution.h
View File

@ -49,6 +49,8 @@ class cxxSolution:public cxxNumKeyword
void Set_cb(LDBLE l_cb) {this->cb = l_cb;}
LDBLE Get_density() const {return this->density;}
void Set_density(LDBLE l_density) {this->density = l_density;}
LDBLE Get_viscosity() const { return this->viscosity; }
void Set_viscosity(LDBLE l_viscos) { this->viscosity = l_viscos; }
LDBLE Get_mass_water() const {return this->mass_water;}
void Set_mass_water(LDBLE l_mass_water) {this->mass_water = l_mass_water;}
LDBLE Get_total_alkalinity() const {return this->total_alkalinity;}
@ -131,6 +133,7 @@ class cxxSolution:public cxxNumKeyword
LDBLE cb;
LDBLE mass_water;
LDBLE density;
LDBLE viscosity;
LDBLE soln_vol;
LDBLE total_alkalinity;
cxxNameDouble totals;

3
SolutionIsotope.cxx
View File

@ -70,7 +70,8 @@ cxxSolutionIsotope::dump_xml(std::ostream & s_oss, unsigned int indent) const
#ifdef NPP
if (!isnan(this->ratio_uncertainty))
#else
if (this->ratio_uncertainty != NAN)
//if (this->ratio_uncertainty != NAN)
if (!std::isnan(this->ratio_uncertainty))
#endif
{
s_oss << indent1;

31
Surface.cxx
View File

@ -36,6 +36,7 @@ cxxSurface::cxxSurface(PHRQ_io *io)
dl_type = NO_DL;
sites_units = SITES_ABSOLUTE;
only_counter_ions = false;
correct_GC = false;
thickness = 1e-8;
debye_lengths = 0.0;
DDL_viscosity = 1.0;
@ -55,6 +56,7 @@ cxxNumKeyword(io)
dl_type = NO_DL;
sites_units = SITES_ABSOLUTE;
only_counter_ions = false;
correct_GC = false;
thickness = 1e-8;
debye_lengths = 0.0;
DDL_viscosity = 1.0;
@ -128,6 +130,8 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons
s_oss << indent1;
s_oss << "-only_counter_ions " << this->only_counter_ions << "\n";
s_oss << indent1;
s_oss << "-correct_GC " << this->correct_GC << "\n";
s_oss << indent1;
s_oss << "-thickness " << this->thickness << "\n";
s_oss << indent1;
s_oss << "-debye_lengths " << this->debye_lengths << "\n";
@ -189,6 +193,7 @@ cxxSurface::read_raw(CParser & parser, bool check)
this->Set_tidied(true);
bool only_counter_ions_defined(false);
//bool correct_GC_defined(false);
bool thickness_defined(false);
bool type_defined(false);
bool dl_type_defined(false);
@ -468,6 +473,17 @@ cxxSurface::read_raw(CParser & parser, bool check)
PHRQ_io::OT_CONTINUE);
}
break;
case 19: // correct_GC
if (!(parser.get_iss() >> this->correct_GC))
{
this->correct_GC = false;
parser.incr_input_error();
parser.
error_msg("Expected boolean value for correct_GC.",
PHRQ_io::OT_CONTINUE);
}
//correct_GC_defined = true;
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
@ -482,6 +498,13 @@ cxxSurface::read_raw(CParser & parser, bool check)
error_msg("Only_counter_ions not defined for SURFACE_RAW input.",
PHRQ_io::OT_CONTINUE);
}
//if (correct_GC_defined == false)
//{
// parser.incr_input_error();
// parser.
// error_msg("correct_GC not defined for SURFACE_RAW input.",
// PHRQ_io::OT_CONTINUE);
//}
if (thickness_defined == false)
{
parser.incr_input_error();
@ -559,6 +582,7 @@ cxxSurface::add(const cxxSurface & addee_in, LDBLE extensive)
if (this->surface_comps.size() == 0)
{
this->only_counter_ions = addee.only_counter_ions;
this->correct_GC = addee.correct_GC;
this->dl_type = addee.dl_type;
this->type = addee.type;
this->sites_units = addee.sites_units;
@ -730,6 +754,7 @@ cxxSurface::Serialize(Dictionary & dictionary, std::vector < int >&ints,
doubles.push_back(this->debye_lengths);
doubles.push_back(this->DDL_viscosity);
doubles.push_back(this->DDL_limit);
ints.push_back(this->correct_GC ? 1 : 0);
ints.push_back(this->transport ? 1 : 0);
this->totals.Serialize(dictionary, ints, doubles);
ints.push_back(this->solution_equilibria ? 1 : 0);
@ -776,6 +801,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->debye_lengths = doubles[dd++];
this->DDL_viscosity = doubles[dd++];
this->DDL_limit = doubles[dd++];
this->correct_GC = (ints[ii++] != 0);
this->transport = (ints[ii++] != 0);
this->totals.Deserialize(dictionary, ints, doubles, ii, dd);
this->solution_equilibria = (ints[ii++] != 0);
@ -803,6 +829,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
std::vector< std::string >::value_type("solution_equilibria"), // 15
std::vector< std::string >::value_type("n_solution"), // 16
std::vector< std::string >::value_type("totals"), // 17
std::vector< std::string >::value_type("tidied") // 18
};
std::vector< std::string >::value_type("tidied"), // 18
std::vector< std::string >::value_type("correct_gc") // 19
};
const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

4
Surface.h
View File

@ -67,6 +67,9 @@ public:
void Set_DDL_viscosity(LDBLE t) {DDL_viscosity = t;}
LDBLE Get_DDL_limit(void) const {return DDL_limit;}
void Set_DDL_limit(LDBLE t) {DDL_limit = t;}
bool Get_correct_GC(void) const { return correct_GC; }
void Set_correct_GC(bool tf) { correct_GC = tf; }
std::vector<LDBLE> Donnan_factors;
bool Get_transport(void) const {return transport;}
void Set_transport(bool tf) {transport = tf;}
cxxNameDouble & Get_totals() {return this->totals;}
@ -91,6 +94,7 @@ protected:
LDBLE debye_lengths;
LDBLE DDL_viscosity;
LDBLE DDL_limit;
bool correct_GC;
bool transport;
cxxNameDouble totals;
bool solution_equilibria;

215
basicsubs.cpp
View File

@ -187,21 +187,60 @@ diff_c(const char *species_name)
/* ---------------------------------------------------------------------- */
{
class species *s_ptr;
LDBLE g;
LDBLE ka, l_z, Dw, ff, sqrt_mu;
sqrt_mu = sqrt(mu_x);
s_ptr = s_search(species_name);
if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/)
//LDBLE g;
//if (s_ptr != NULL /*&& s_ptr->in != FALSE && s_ptr->type < EMINUS*/)
//{
// g = s_ptr->dw;
// if (s_ptr->dw_t)
// g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
// g *= viscos_0_25 / viscos * tk_x / 298.15;
//}
//else
//{
// g = 0;
//}
//return (g);
if (s_ptr == NULL)
return(0);
if ((Dw = s_ptr->dw) == 0)
{
g = s_ptr->dw;
if (s_ptr->dw_t)
g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
g *= viscos_0_25 / viscos * tk_x / 298.15;
if (correct_Dw)
Dw = default_Dw;
else
return(0);
}
if ((l_z = fabs(s_ptr->z)) == 0)
{
//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
}
else
{
g = 0;
if (s_ptr->dw_a2)
ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
else
ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
if (s_ptr->dw_a)
{
ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
//if (print_viscosity && s_ptr->dw_a_visc)
// ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc);
}
else
{
ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
}
Dw *= ff;
}
return (g);
if (tk_x != 298.15 && s_ptr->dw_t)
Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
s_ptr->dw_corr = Dw;
return (Dw * viscos_0_25 / viscos_0);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
@ -209,22 +248,58 @@ setdiff_c(const char *species_name, double d)
/* ---------------------------------------------------------------------- */
{
class species *s_ptr;
LDBLE g;
LDBLE ka, l_z, Dw, ff, sqrt_mu;
sqrt_mu = sqrt(mu_x);
s_ptr = s_search(species_name);
if (s_ptr != NULL)
//LDBLE g;
//s_ptr = s_search(species_name);
//if (s_ptr != NULL)
//{
// s_ptr->dw = d;
// g = s_ptr->dw;
// if (s_ptr->dw_t)
// g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
// g *= viscos_0_25 / viscos * tk_x / 298.15;;
//}
//else
//{
// g = 0;
//}
//return (g);
if (s_ptr == NULL)
return(0);
Dw = s_ptr->dw = d;
if ((l_z = fabs(s_ptr->z)) == 0)
{
s_ptr->dw = d;
g = s_ptr->dw;
if (s_ptr->dw_t)
g *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
g *= viscos_0_25 / viscos * tk_x / 298.15;;
//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
}
else
{
g = 0;
if (s_ptr->dw_a2)
ka = DH_B * s_ptr->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
else
ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
if (s_ptr->dw_a)
{
ff = exp(-s_ptr->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
//if (print_viscosity && s_ptr->dw_a_visc)
// ff *= pow((viscos_0 / viscos), s_ptr->dw_a_visc);
}
else
{
ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
}
Dw *= ff;
}
return (g);
if (tk_x != 298.15 && s_ptr->dw_t)
Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
s_ptr->dw_corr = Dw;
return (Dw * viscos_0_25 / viscos_0);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
@ -314,19 +389,18 @@ calc_SC(void)
SC = 0;
//LDBLE ta1, ta2, ta3, ta4;
for (i = 0; i < (int)this->s_x.size(); i++)
{
// ** for optimizing, get numbers from -analyt for H+ = H+...
//if (!strcmp(s_x[i]->name, "H+"))
//{
// ta1 = s_x[i]->logk[2];
// ta2 = s_x[i]->logk[3];
// ta3 = s_x[i]->logk[4];
// ta4 = s_x[i]->logk[5];
// break;
//}
//
}
//for (i = 0; i < (int)this->s_x.size(); i++)
//{
// // ** for optimizing, get numbers from -analyt for H+ = H+...
// if (!strcmp(s_x[i]->name, "H+"))
// {
// ta1 = s_x[i]->logk[2] * 1e15;
// ta2 = s_x[i]->logk[3] * 1e15;
// ta3 = s_x[i]->logk[4] * 1e15;
// ta4 = s_x[i]->logk[5] * 1e15;
// break;
// }
//}
for (i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type != AQ && s_x[i]->type != HPLUS)
@ -338,43 +412,51 @@ calc_SC(void)
else
continue;
}
if (s_x[i]->lm < min_dif_LM)
continue;
if ((l_z = fabs(s_x[i]->z)) == 0)
l_z = 1; // only a 1st approximation for correct_Dw in electrical field
if (s_x[i]->dw_t)
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); // the viscosity multiplier is done in SC
if (s_x[i]->dw_a2)
ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
{
//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
}
else
{
ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x , 0.75));
//ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2));
//ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2);
}
if (s_x[i]->dw_a)
{
ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
if (print_viscosity)
if (s_x[i]->dw_a2)
ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
else
{
ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc);
ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.75));
//ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2));
//ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2);
}
if (s_x[i]->dw_a)
{
ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
//if (print_viscosity && s_x[i]->dw_a_visc)
// ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc);
}
else
{
ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
//ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka));
}
Dw *= ff;
}
else
if (tk_x != 298.15)
{
ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
//ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka));
if (s_x[i]->dw_t)
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15);
//else
//{
// Dw *= exp(ta1 / tk_x - ta1 / 298.15);
//}
}
Dw *= ff;
s_x[i]->dw_corr = Dw;
if (s_x[i]->z == 0)
continue;
s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw;
SC += s_x[i]->dw_t_SC;
}
SC *= 1e7 * F_C_MOL * F_C_MOL / (R_KJ_DEG_MOL * 298150.0);
/* correct for temperature dependency...
SC_T = SC_298 * (Dw_T / T) * (298 / Dw_298) and
Dw_T = Dw_298 * (T / 298) * (viscos_298 / viscos_T) give:
/* correct for viscosity dependency...
SC_T = SC_298 * (viscos_298 / viscos_T)
*/
SC *= viscos_0_25 / viscos_0;
@ -1520,7 +1602,8 @@ get_calculate_value(const char *name)
#ifdef NPP
if (isnan(rate_moles))
#else
if (rate_moles == NAN)
//if (rate_moles == NAN)
if(std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
@ -1721,22 +1804,22 @@ pressure(void)
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
pr_phi(const char* phase_name)
pr_phi(const char *phase_name)
/* ---------------------------------------------------------------------- */
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
if (gas_phase_ptr != NULL)
{
int l;
int l;
class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE);
if (phase_ptr == NULL)
{
error_string = sformatf("Gas %s, not found.", phase_name);
warning_msg(error_string);
return (1e-99);
}
if (phase_ptr == NULL)
{
error_string = sformatf( "Gas %s, not found.", phase_name);
warning_msg(error_string);
return (1e-99);
}
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]);
int j;
class phase* phase_ptr_gas = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
@ -1744,8 +1827,8 @@ pr_phi(const char* phase_name)
{
if (gas_phase_ptr->Get_pr_in())
{
return phase_ptr->pr_phi;
}
return phase_ptr->pr_phi;
}
else
{
return gas_comp_ptr->Get_phi();
@ -1753,7 +1836,7 @@ pr_phi(const char* phase_name)
}
}
}
return(1.0);
return (1.0);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::

3
cl1.cpp
View File

@ -1,6 +1,7 @@
#include <stdio.h>
#include <string.h>
#include <math.h>
//#include <math.h>
#include <cmath>
#include <stdlib.h>
#include "Phreeqc.h"
#include "phqalloc.h"

4
class_main.cpp
View File

@ -323,7 +323,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = snprintf(buffer, sizeof(buffer), "%s", "July 5, 2021");
len = snprintf(buffer, sizeof(buffer), "%s", "March 20, 2023");
#else
len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@");
#endif
@ -507,7 +507,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.7.1, compiled July 5, 2021\n"));
output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 20, 2023\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token.c_str()));
#ifdef NPP

3
common/Utils.cxx
View File

@ -10,7 +10,8 @@
#include "Utils.h"
#include "Parser.h"
#include "float.h"
#include "math.h"
//#include <math.h>
#include <cmath>
#if defined(PHREEQCI_GUI)
#ifdef _DEBUG

2
gases.cpp
View File

@ -624,6 +624,7 @@ int Phreeqc::
calc_fixed_volume_gas_pressures(void)
/* ---------------------------------------------------------------------- */
{
//int n_g = 0;
LDBLE lp;
class rxn_token *rxn_ptr;
class phase *phase_ptr;
@ -644,6 +645,7 @@ calc_fixed_volume_gas_pressures(void)
{
if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0)
PR = true;
//n_g++;
}
gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gas_unknowns[i]->moles);
}

4
global_structures.h
View File

@ -5,8 +5,8 @@
/* ----------------------------------------------------------------------
* #define DEFINITIONS
* ---------------------------------------------------------------------- */
#ifndef NAN
# define NAN -99999999
#if !defined(NAN)
#define NAN nan("1")
#endif
#define MISSING -9999.999
#include "NA.h" /* NA = not available */

407
integrate.cpp
View File

@ -25,9 +25,9 @@ calc_all_g(void)
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* calculate g for each surface
*/
/*
* calculate g for each surface
*/
epsilon = convergence_tolerance;
if (convergence_tolerance >= 1e-8)
{
@ -45,7 +45,7 @@ calc_all_g(void)
if (x[j]->type != SURFACE_CB)
continue;
if (debug_diffuse_layer == TRUE)
output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
output_msg(sformatf("Calc_all_g, X[%d]\n", j));
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
std::map<LDBLE, cxxSurfDL> temp_g_map;
cxxSurfDL temp_g;
@ -56,11 +56,11 @@ calc_all_g(void)
/* 1000 J/kJ and 1000 L/m**3 */
//alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
// tk_x * 0.5);
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
tk_x * 0.5);
/*
* calculate g for given surface for each species
*/
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
tk_x * 0.5);
/*
* calculate g for given surface for each species
*/
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
@ -70,10 +70,10 @@ calc_all_g(void)
z_global = s_x[i]->z;
if (charge_ptr->Get_grams() > 0.0)
{
if ((use.Get_surface_ptr()->Get_only_counter_ions() == false) ||
(((x[j]->master[0]->s->la > 0) && (z_global < 0))
|| ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
|| ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
{
if (xd_global > 0.1)
{
@ -180,11 +180,11 @@ calc_all_g(void)
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double) z_global,
(double) charge_ptr->Get_g_map()[z_global].Get_g(),
(double) new_g,
(double) (new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double)z_global,
(double)charge_ptr->Get_g_map()[z_global].Get_g(),
(double)new_g,
(double)(new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
}
}
charge_ptr->Get_g_map()[z_global].Set_g(new_g);
@ -201,7 +201,7 @@ calc_all_g(void)
g_function(xd_global) / F_C_MOL;
dg *=
-2. / (exp(x[j]->master[0]->s->la * LOG_10) *
exp(x[j]->master[0]->s->la * LOG_10));
exp(x[j]->master[0]->s->la * LOG_10));
if ((xd_global - 1) < 0.0)
{
dg *= -1.0;
@ -225,16 +225,16 @@ calc_all_g(void)
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf("\nSurface component %d: charge,\tg,\tdg/dlny,\txd\n",
(int) charge_ptr->Get_g_map().size()));
(int)charge_ptr->Get_g_map().size()));
std::map<LDBLE, cxxSurfDL>::iterator it;
for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
{
output_msg(sformatf(
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double) it->first,
(double) it->second.Get_g(),
(double) it->second.Get_dg(),
(double) xd_global));
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double)it->first,
(double)it->second.Get_g(),
(double)it->second.Get_dg(),
(double)xd_global));
}
}
}
@ -253,10 +253,10 @@ g_function(LDBLE x_value)
return (0.0);
sum = 0.0;
ln_x_value = log(x_value);
cxxSurfaceCharge *charge_ptr = &(use.Get_surface_ptr()->Get_surface_charges()[0]);
std::map<LDBLE, cxxSurfDL>::iterator it = charge_ptr->Get_g_map().begin();
for ( ; it != charge_ptr->Get_g_map().end(); it++)
for (; it != charge_ptr->Get_g_map().end(); it++)
{
it->second.Set_psi_to_z(exp(ln_x_value * it->first) - 1.0);
}
@ -272,31 +272,31 @@ g_function(LDBLE x_value)
sum = 0.0;
sum1 = 0.0;
output_msg(sformatf(
"Species\tmoles\tX**z-1\tsum\tsum charge\n"));
"Species\tmoles\tX**z-1\tsum\tsum charge\n"));
for (i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
{
sum += s_x[i]->moles * (pow(x_value, s_x[i]->z) - 1.0);
sum1 += s_x[i]->moles * s_x[i]->z;
output_msg(sformatf( "%s\t%e\t%e\t%e\t%e\n",
s_x[i]->name, (double) s_x[i]->moles,
(double) (pow((LDBLE) x_value, (LDBLE) s_x[i]->z) -
1.0), (double) sum, (double) sum1));
output_msg(sformatf("%s\t%e\t%e\t%e\t%e\n",
s_x[i]->name, (double)s_x[i]->moles,
(double)(pow((LDBLE)x_value, (LDBLE)s_x[i]->z) -
1.0), (double)sum, (double)sum1));
}
}
error_string = sformatf( "Negative sum in g_function, %e\t%e.",
(double) sum, (double) x_value);
error_string = sformatf("Negative sum in g_function, %e\t%e.",
(double)sum, (double)x_value);
error_msg(error_string, CONTINUE);
error_string = sformatf(
"Solutions must be charge balanced, charge imbalance is %e\n",
(double) sum1);
"Solutions must be charge balanced, charge imbalance is %e\n",
(double)sum1);
error_msg(error_string, STOP);
}
return_value =
(exp(ln_x_value * z_global) -
1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
return (return_value);
}
/* ---------------------------------------------------------------------- */
@ -309,9 +309,9 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
ns = 1;
dif = fabs(xv - xa[1]);
/*
* Malloc work space
*/
/*
* Malloc work space
*/
std::vector<double> c, d;
c.resize((size_t)n + 1);
d.resize((size_t)n + 1);
@ -354,7 +354,7 @@ polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
}
*yv += *dy;
/* *yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
/* *yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
}
return;
}
@ -376,7 +376,7 @@ midpnt(LDBLE x1, LDBLE x2, int n)
{
for (it = 1, j = 1; j < n - 1; j++)
it *= 3;
tnm = (LDBLE) it;
tnm = (LDBLE)it;
del = (x2 - x1) / (3 * tnm);
ddel = del + del;
xv = x1 + 0.5 * del;
@ -420,7 +420,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"Iterations in qromb_midpnt: %d\n", j));
"Iterations in qromb_midpnt: %d\n", j));
}
return (sv[j]);
}
@ -436,7 +436,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"Iterations in qromb_midpnt: %d\n", j));
"Iterations in qromb_midpnt: %d\n", j));
}
return (ss);
}
@ -444,7 +444,7 @@ qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
}
error_string = sformatf(
"\nToo many iterations integrating diffuse layer.\n");
"\nToo many iterations integrating diffuse layer.\n");
error_msg(error_string, STOP);
return (-999.9);
}
@ -455,9 +455,9 @@ calc_init_g(void)
{
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* calculate g for each surface
*/
/*
* calculate g for each surface
*/
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
@ -468,7 +468,7 @@ calc_init_g(void)
/* second 1000 is liters/m**3 */
//alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// 1000.0 * tk_x * 0.5);
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
1000.0 * tk_x * 0.5);
if (charge_ptr->Get_g_map().size() == 0)
@ -476,9 +476,9 @@ calc_init_g(void)
cxxSurfDL temp_g;
charge_ptr->Get_g_map()[0.0] = temp_g;
}
/*
* calculate g for given surface for each species
*/
/*
* calculate g for given surface for each species
*/
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
@ -511,7 +511,7 @@ calc_init_g(void)
{
int is = s_x[i]->number;
assert (is < (int) s_diff_layer.size());
assert(is < (int)s_diff_layer.size());
// species found in diff_layer
s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0);
s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0);
@ -520,15 +520,15 @@ calc_init_g(void)
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"\nSurface component %d: charge,\tg,\tdg\n",
(int) charge_ptr->Get_g_map().size()));
"\nSurface component %d: charge,\tg,\tdg\n",
(int)charge_ptr->Get_g_map().size()));
std::map<LDBLE, cxxSurfDL>::iterator it;
for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) it->first,
(double) it->second.Get_g(),
(double) it->second.Get_dg()));
output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\n",
(double)it->first,
(double)it->second.Get_g(),
(double)it->second.Get_dg()));
}
}
}
@ -539,21 +539,21 @@ int Phreeqc::
initial_surface_water(void)
/* ---------------------------------------------------------------------- */
{
/*
* In initial surface calculation, need to calculate
* mass of water in diffuse layer.
* diffuse layer water + aqueous solution water = bulk water.
* Ionic strength is fixed, so diffuse-layer water will not change
*/
/*
* In initial surface calculation, need to calculate
* mass of water in diffuse layer.
* diffuse layer water + aqueous solution water = bulk water.
* Ionic strength is fixed, so diffuse-layer water will not change
*/
LDBLE debye_length, b, r, rd, ddl_limit, rd_limit, fraction, sum_surfs, l_s;
LDBLE damp_aq;
/*
* Debye length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
*
* 1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
*/
//debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
// / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
/*
* Debye length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
*
* 1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
*/
//debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
// / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
debye_length = (eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
debye_length = sqrt(debye_length);
@ -561,10 +561,10 @@ initial_surface_water(void)
/* ddl is at most the fraction ddl_limit of bulk water */
ddl_limit = use.Get_surface_ptr()->Get_DDL_limit();
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
if (use.Get_surface_ptr()->Get_debye_lengths() > 0)
{
@ -606,14 +606,14 @@ initial_surface_water(void)
mass_water_surfaces_x =
use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 - ddl_limit);
r = 0.002 * (mass_water_surfaces_x +
use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
rd = rd_limit;
use.Get_surface_ptr()->Set_thickness(rd);
}
else
mass_water_surfaces_x =
(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
for (int i = 0; i < count_unknowns; i++)
{
if (x[i]->type != SURFACE_CB)
@ -692,14 +692,14 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* Find position of component in list of components
*/
/*
* Find position of component in list of components
*/
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
count_elts = 0;
paren_count = 0;
mass_water_surface = charge_ptr->Get_mass_water();
@ -717,9 +717,9 @@ sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
}
moles_excess = mass_water_aq_x * molality * g;
moles_surface = mass_water_surface * molality + moles_excess;
/*
* Accumulate elements in diffuse layer
*/
/*
* Accumulate elements in diffuse layer
*/
add_elt_list(s_x[j]->next_elt, moles_surface);
}
add_elt_list(s_h2o->next_elt, mass_water_surface / gfw_water);
@ -731,22 +731,33 @@ int Phreeqc::
calc_all_donnan(void)
/* ---------------------------------------------------------------------- */
{
bool converge;
bool converge;
int cd_m;
LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot;
LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot, co_ion;
LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1;
if (use.Get_surface_ptr() == NULL)
return (OK);
//f_sinh = sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// tk_x * mu_x);
f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
tk_x * mu_x);
/*
* calculate g for each surface...
*/
if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths()) // DL_pitz
tk_x * mu_x);
bool only_count = use.Get_surface_ptr()->Get_only_counter_ions();
bool correct_GC = use.Get_surface_ptr()->Get_correct_GC();
/* calculate g for each surface...
*/
if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths() ||
correct_GC) // DL_pitz && correct_GC
initial_surface_water();
LDBLE nDbl = 1;
if (correct_GC)
{
if ((nDbl = use.Get_surface_ptr()->Get_debye_lengths()) == 0)
{
LDBLE debye_length = f_sinh / (F_C_MOL * mu_x * 4e3);
nDbl = use.Get_surface_ptr()->Get_thickness() / debye_length;
//use.Get_surface_ptr()->Set_debye_lengths(nDbl);
}
}
converge = TRUE;
for (int j = 0; j < count_unknowns; j++)
{
@ -755,10 +766,10 @@ calc_all_donnan(void)
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
/*
* sum eq of each charge number in solution...
*/
output_msg(sformatf("Calc_all_g, X[%d]\n", j));
/*
* sum eq of each charge number in solution...
*/
std::map<LDBLE, LDBLE>::iterator it;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
@ -778,15 +789,18 @@ calc_all_donnan(void)
f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2;
cd_m = 1;
} else
}
else
{
f_psi = x[j]->master[0]->s->la * LOG_10;
cd_m = -1;
}
surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
if (fabs(surf_chrg_eq) > 5e3)
LDBLE lim_seq = 5e3;
if (correct_GC) lim_seq = 5e1;
if (fabs(surf_chrg_eq) > lim_seq)
{
surf_chrg_eq = (surf_chrg_eq < 0 ? -5e3 : 5e3);
surf_chrg_eq = (surf_chrg_eq < 0 ? -lim_seq : lim_seq);
var1 = surf_chrg_eq / (A_surf * f_sinh / F_C_MOL);
var1 = (var1 + sqrt(var1 * var1 + 1));
f_psi = (var1 > 1e-8 ? log(var1) : -18.4);
@ -799,8 +813,11 @@ calc_all_donnan(void)
surf_chrg_eq2 = A_surf * f_sinh * sinh(f_psi2) / F_C_MOL;
/* find psi_avg that matches surface charge... */
psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq);
psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2);
std::vector<LDBLE> zcorr(charge_group_map.size());
std::vector<LDBLE> zcorr2(charge_group_map.size());
//LDBLE fD = 0;
psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq, nDbl, zcorr);
psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2, nDbl, zcorr2);
/*output_msg(sformatf( "psi's %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */
@ -808,9 +825,15 @@ calc_all_donnan(void)
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
ratio_aq_tot = charge_ptr->Get_mass_water() / mass_water_bulk_x;
int z_iter = 0;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
LDBLE z = it->first;
LDBLE z = it->first, z1 = z;
co_ion = surf_chrg_eq * z;
if (correct_GC)
z1 = zcorr[z_iter];
//z1 *= cgc[0] * pow(z_factor, abs(z));
if (!ratio_aq)
{
charge_ptr->Get_g_map()[z].Set_g(0);
@ -819,17 +842,13 @@ calc_all_donnan(void)
converge = true;
continue;
}
new_g = ratio_aq * (exp(cd_m * z * psi_avg) - 1);
if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
//((surf_chrg_eq < 0 && z < 0)
// || (surf_chrg_eq > 0 && z > 0)))
new_g = ratio_aq * (exp(cd_m * z1 * psi_avg) - 1);
if (only_count && co_ion > 0)
new_g = -ratio_aq;
if (new_g <= -ratio_aq)
new_g = -ratio_aq + G_TOL * 1e-3;
new_g2 = ratio_aq * (exp(cd_m * z * psi_avg2) - 1);
if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
//((surf_chrg_eq < 0 && z < 0)
// || (surf_chrg_eq > 0 && z > 0)))
new_g2 = ratio_aq * (exp(cd_m * z1 * psi_avg2) - 1);
if (only_count && co_ion > 0)
new_g2 = -ratio_aq;
if (new_g2 <= -ratio_aq)
new_g2 = -ratio_aq + G_TOL * 1e-3;
@ -859,9 +878,9 @@ calc_all_donnan(void)
/* save Boltzmann factor * water fraction for MCD calc's in transport */
if (converge)
{
if (use.Get_surface_ptr()->Get_only_counter_ions())
if (only_count)
{
if (surf_chrg_eq * z > 0) // co-ions are not in the DL
if (co_ion > 0) // co-ions are not in the DL
charge_ptr->Get_z_gMCD_map()[z] = 0;
else // assume that counter-ions have the free water conc for diffusion
charge_ptr->Get_z_gMCD_map()[z] = ratio_aq_tot;
@ -869,23 +888,27 @@ calc_all_donnan(void)
else
charge_ptr->Get_z_gMCD_map()[z] = (new_g / ratio_aq + 1) * ratio_aq_tot;
}
z_iter++;
}
if (debug_diffuse_layer == TRUE)
{
std::string name = x[j]->master[0]->elt->name;
std::string name = x[j]->master[0]->elt->name;
Utilities::replace("_psi", "", name);
output_msg(sformatf(
"\nDonnan all on %s (%d): charge, \tg, \tdg, Psi_surface = %8f V. \n",
name.c_str(), (int) charge_ptr->Get_g_map().size(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
tk_x / F_KJ_V_EQ));
"\nDonnan all on %s (%d): charge, \tg, \tdg, \tzcorr, \tPsi_surface = %8f V, \tDebye lengths = %8f. \n",
name.c_str(), (int)charge_ptr->Get_g_map().size(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ,
nDbl));
int i = 0;
for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
i_it != charge_ptr->Get_g_map().end(); i_it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) i_it->first,
(double) i_it->second.Get_g(),
(double) i_it->second.Get_dg()));
output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double)i_it->first,
(double)i_it->second.Get_g(),
(double)i_it->second.Get_dg(),
(double)zcorr[i]));
i++;
}
}
}
@ -904,7 +927,7 @@ calc_init_donnan(void)
//sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// tk_x * mu_x);
sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
tk_x * mu_x);
tk_x * mu_x);
if (convergence_tolerance >= 1e-8)
{
G_TOL = 1e-9;
@ -913,9 +936,9 @@ calc_init_donnan(void)
{
G_TOL = 1e-13;
}
/*
* sum eq of each charge number in solution...
*/
/*
* sum eq of each charge number in solution...
*/
charge_group_map.clear();
charge_group_map[0.0] = 0.0;
@ -932,9 +955,10 @@ calc_init_donnan(void)
charge_group_map[s_x[i]->z] = s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
}
}
/*
* calculate g for each surface...
*/
std::vector<LDBLE> zcorr(charge_group_map.size());
/*
* calculate g for each surface...
*/
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
@ -948,14 +972,13 @@ calc_init_donnan(void)
{
f_psi = x[(size_t)j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2;
} else
}
else
f_psi = x[j]->master[0]->s->la * LOG_10;
surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
/* find psi_avg that matches surface charge... */
/* psi_avg = calc_psi_avg (0);
appt 7/9/8... may have to change above one */
psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq);
psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq, 0, zcorr);
/* fill in g's */
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
@ -975,7 +998,7 @@ calc_init_donnan(void)
if (charge_ptr->Get_g_map()[z].Get_g() != 0)
{
charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) /
charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) /
(eq * F_C_MOL));
}
else
@ -986,7 +1009,7 @@ calc_init_donnan(void)
for (int i = 0; i < (int)this->s_x.size(); i++)
{
int is = s_x[i]->number;
assert (is < (int) s_diff_layer.size());
assert(is < (int)s_diff_layer.size());
s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0.0);
s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0.0);
@ -997,17 +1020,17 @@ calc_init_donnan(void)
std::string name = x[j]->master[0]->elt->name;
Utilities::replace("_psi", "", name);
output_msg(sformatf(
"\nDonnan init on %s : charge, \tg, \tdg, Psi_surface = %8f V. \n",
name.c_str(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
tk_x / F_KJ_V_EQ));
"\nDonnan init on %s : charge, \tg, \tdg, \tzcorr, Psi_surface = %8f V. \n",
name.c_str(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
tk_x / F_KJ_V_EQ));
for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
i_it != charge_ptr->Get_g_map().end(); i_it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) i_it->first,
(double) i_it->second.Get_g(),
(double) i_it->second.Get_dg()));
output_msg(sformatf("\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double)i_it->first,
(double)i_it->second.Get_g(),
(double)i_it->second.Get_dg()));
}
}
}
@ -1015,13 +1038,13 @@ calc_init_donnan(void)
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq, LDBLE nDbl, std::vector<LDBLE> &zcorr)
/* ---------------------------------------------------------------------- */
{
/*
* calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
*/
LDBLE fd, fd1, p, temp, ratio_aq;
/*
* calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
*/
LDBLE fd, fd1, p, /*psi_DL, */p_psi = R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ, temp, ratio_aq, z, z1, z1_c, eq, co_ion, sum_counter, sum_co;
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
p = 0;
@ -1031,30 +1054,79 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
p = -0.5 * log(-surf_chrg_eq * ratio_aq / mu_x + 1);
else if (surf_chrg_eq > 0)
p = 0.5 * log(surf_chrg_eq * ratio_aq / mu_x + 1);
/*
* Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
* g(p) = exp(-p * z_i) * ratio_aq
* Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
* g(p) = (exp(-p *z_i) - 1) * ratio_aq
*/
int l_iter = 0;
/*
* Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
* g(p) = exp(-p * z_i) * ratio_aq
* Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
* g(p) = (exp(-p *z_i) - 1) * ratio_aq
* with correct_GC true:
* correct ions to better match the integrated concentrations:
z == 1? z *= 0.285 cgc[6]
z == 2? z *= 0.372 cgc[7]
z == -1? z *= cgc[0] * (mu_x**( cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4])
z == -2? z *= cgc[0] * (mu_x**(cgc[5] * cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * I ** cgc[4])
*/
cxxSurface *surf_p = use.Get_surface_ptr();
bool correct_GC = surf_p->Get_correct_GC(), local_correct_GC = correct_GC;
bool only_count = surf_p->Get_only_counter_ions();
LDBLE Gamma = fabs(surf_chrg_eq) / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()) / 1e-6,
cgc[10] = { 0.36, 0.1721, 0.798, 0.287, 0.1457, 1.2, 0.285, 0.372 };
if (!surf_p->Donnan_factors.empty())
std::copy(std::begin(surf_p->Donnan_factors), std::end(surf_p->Donnan_factors), cgc);
cgc[1] *= pow(nDbl, cgc[2]) * pow(Gamma, cgc[3]) * pow(mu_x, cgc[4]);
int l_iter = 0, z_iter;
sum_co = sum_counter = 0;
do
{
fd = surf_chrg_eq;
fd1 = 0.0;
z_iter = 0;
if (l_iter == 1 && local_correct_GC && fabs(sum_counter) < fabs(sum_co))
{
local_correct_GC = false;
l_iter = 0;
}
std::map<LDBLE, LDBLE>::iterator it;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
LDBLE z = it->first;
if (!z || (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0))
z = it->first;
z1 = z;
if (l_iter == 0) zcorr[z_iter] = z;
co_ion = surf_chrg_eq * z;
if (!z || (only_count && co_ion > 0))
{
z_iter++;
continue;
LDBLE eq = it->second;
/* multiply with ratio_aq for multiplier options cp and cm
in calc_all_donnan (not used now)... */
temp = exp(-z * p) * ratio_aq;
}
if (nDbl && local_correct_GC)
{
/*psi_DL = fabs(p * p_psi);*/
if (co_ion < 0)
{//counter-ion
if (fabs(z) > 1) temp = cgc[7];
else temp = cgc[6];
sum_counter += z * temp;
}
else
{// co-ion
if (fabs(z) > 1) temp = cgc[0] * pow(mu_x, cgc[1] * cgc[5]);
else temp = cgc[0] * pow(mu_x, cgc[1]);
sum_co += z * temp;
}
zcorr[z_iter] = z * temp;
}
z1 = zcorr[z_iter];
eq = it->second;
temp = exp(-z1 * p) * ratio_aq;
fd += eq * temp;
fd1 -= z * eq * temp;
fd1 -= z1 * eq * temp;
if (z == 1) z1_c = z1;
z_iter++;
}
fd /= -fd1;
p += (fd > 1) ? 1 : ((fd < -1) ? -1 : fd);
@ -1076,12 +1148,11 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
error_msg(error_string, STOP);
}
}
while (fabs(fd) > 1e-12 && p != 0.0);
} while (fabs(fd) > 1e-12 && p != 0.0);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf(
"iter in calc_psi_avg = %d. g(+1) = %8f. surface charge = %12.4e.\n",
l_iter, (double) (exp(-p) - 1), (double) surf_chrg_eq));
"iter in calc_psi_avg = %d. g(+1) = %8f, surface charge = %12.4e, psi_DL = %12.3e V.\n",
l_iter, (double)(exp(-p) - 1), (double)surf_chrg_eq, (double)(p * z1_c * p_psi)));
return (p);
}

14
inverse.cpp
View File

@ -3586,8 +3586,10 @@ check_isotopes(class inverse *inv_ptr)
if (j < inv_ptr->i_u[i].uncertainties.size()
&& !isnan(inv_ptr->i_u[i].uncertainties[j]))
#else
//if (j < inv_ptr->i_u[i].uncertainties.size()
// && inv_ptr->i_u[i].uncertainties[j] != NAN)
if (j < inv_ptr->i_u[i].uncertainties.size()
&& inv_ptr->i_u[i].uncertainties[j] != NAN)
&& !std::isnan(inv_ptr->i_u[i].uncertainties[j]))
#endif
{
kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]);
@ -3598,8 +3600,10 @@ check_isotopes(class inverse *inv_ptr)
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
// && inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN)
else if (inv_ptr->i_u[i].uncertainties.size() > 0
&& inv_ptr->i_u[i].uncertainties[(size_t)inv_ptr->i_u[i].uncertainties.size() - 1] != NAN)
&& !std::isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
#endif
{
kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]);
@ -3609,7 +3613,8 @@ check_isotopes(class inverse *inv_ptr)
#ifdef NPP
else if (!isnan(kit->second.Get_ratio_uncertainty()))
#else
else if (kit->second.Get_ratio_uncertainty() != NAN)
//else if (kit->second.Get_ratio_uncertainty() != NAN)
else if (!std::isnan(kit->second.Get_ratio_uncertainty()))
#endif
{
kit->second.Set_x_ratio_uncertainty(
@ -3641,7 +3646,8 @@ check_isotopes(class inverse *inv_ptr)
#ifdef NPP
if (isnan(kit->second.Get_x_ratio_uncertainty()))
#else
if (kit->second.Get_x_ratio_uncertainty() == NAN)
//if (kit->second.Get_x_ratio_uncertainty() == NAN)
if (std::isnan(kit->second.Get_x_ratio_uncertainty()))
#endif
{
error_string = sformatf(

9
isotopes.cpp
View File

@ -903,7 +903,8 @@ punch_calculate_values(void)
#ifdef NPP
if (isnan(rate_moles))
#else
if (rate_moles == NAN)
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
@ -1136,7 +1137,8 @@ calculate_values(void)
#ifdef NPP
if (isnan(rate_moles))
#else
if (rate_moles == NAN)
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
@ -1203,7 +1205,8 @@ calculate_values(void)
#ifdef NPP
if (isnan(rate_moles))
#else
if (rate_moles == NAN)
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",

3
kinetics.cpp
View File

@ -118,7 +118,8 @@ calc_kinetic_reaction(cxxKinetics *kinetics_ptr, LDBLE time_step)
#ifdef NPP
if (isnan(rate_moles))
#else
if (rate_moles == NAN)
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Moles of reaction not SAVEed for %s.",

10
mainsubs.cpp
View File

@ -540,6 +540,8 @@ initial_exchangers(int print)
viscosity();
species_list_sort();
print_exchange();
if (pr.user_print)
print_user_print();
xexchange_save(n_user);
punch_all();
/* free_model_allocs(); */
@ -671,7 +673,9 @@ initial_gas_phases(int print)
}
print_gas_phase();
if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/)
if (pr.user_print)
print_user_print();
if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/)
warning_msg("While initializing gas phase composition by equilibrating:\n"
" Found definitions of gas` critical temperature and pressure.\n"
" Going to use Peng-Robinson in subsequent calculations.\n");
@ -745,7 +749,8 @@ initial_surfaces(int print)
set_and_run_wrapper(-1, FALSE, FALSE, -1, 0.0);
species_list_sort();
print_surface();
/*print_all(); */
if (pr.user_print)
print_user_print();
punch_all();
xsurface_save(n_user);
/* free_model_allocs(); */
@ -1257,6 +1262,7 @@ xsolution_save(int n_user)
temp_solution.Set_ah2o(ah2o_x);
//temp_solution.Set_density(density_x);
temp_solution.Set_density(calc_dens());
temp_solution.Set_viscosity(this->viscos);
temp_solution.Set_total_h(total_h_x);
temp_solution.Set_total_o(total_o_x);
temp_solution.Set_cb(cb_x); /* cb_x does not include surface charge sfter sum_species */

5
model.cpp
View File

@ -2543,6 +2543,7 @@ int Phreeqc::
calc_gas_pressures(void)
/* ---------------------------------------------------------------------- */
{
//int n_g = 0;
LDBLE lp, V_m = 0;
class rxn_token *rxn_ptr;
std::vector<class phase *> phase_ptrs;
@ -2575,6 +2576,7 @@ calc_gas_pressures(void)
phase_ptrs.push_back(phase_ptr);
if (!PR && phase_ptr->t_c > 0 && phase_ptr->p_c > 0)
PR = true;
//n_g++;
}
if (iterations > 2 && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME)
{
@ -3776,6 +3778,7 @@ residuals(void)
int converge;
LDBLE l_toler;
//LDBLE sum_residual;
LDBLE sinh_constant;
LDBLE sum, sum1;
class master *master_ptr, *master_ptr1, *master_ptr2;
@ -3783,6 +3786,7 @@ residuals(void)
int print_fail;
std::vector<LDBLE> cd_psi;
print_fail = FALSE;
//sum_residual = 0.0;
sigmaddl = 0;
sum = 0;
/*
@ -4526,6 +4530,7 @@ residuals(void)
* Store residuals in array
*/
my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i];
//sum_residual += fabs(residual[i]);
}
/*
* Return

18
pitzer.cpp
View File

@ -52,6 +52,7 @@ pitzer_tidy(void)
const char *string1, *string2;
int i, j, order;
int i0, i1, i2;
//int count_pos, count_neg, count_neut, count[3], jj;
int count_neut, count[3], jj;
LDBLE z0, z1;
class pitz_param *pzp_ptr;
@ -339,21 +340,22 @@ pitzer_tidy(void)
i0 = pitz_params[i]->ispec[0];
i1 = pitz_params[i]->ispec[1];
i2 = pitz_params[i]->ispec[2];
//count_pos = count_neg = count_neut = 0;
count_neut = 0;
for (j = 0; j <= 2; j++)
{
//if (spec[pitz_params[i]->ispec[j]]->z > 0)
//{
// count_pos++;
//}
if (spec[pitz_params[i]->ispec[j]]->z > 0)
{
//count_pos++;
}
if (spec[pitz_params[i]->ispec[j]]->z == 0)
{
count_neut++;
}
//if (spec[pitz_params[i]->ispec[j]]->z < 0)
//{
// count_neg++;
//}
if (spec[pitz_params[i]->ispec[j]]->z < 0)
{
//count_neg++;
}
}
/* All neutral */
if (count_neut == 3)

5
prep.cpp
View File

@ -1980,6 +1980,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr)
* Output: space is allocated and a list of master species pointers is
* returned.
*/
//int j, l, count_list;
int j, l;
char token[MAX_LENGTH];
std::vector<class master*> master_ptr_list;
@ -1987,6 +1988,7 @@ get_list_master_ptrs(const char* cptr, class master *master_ptr)
/*
* Make list of master species pointers
*/
//count_list = 0;
//master_ptr_list = unknown_alloc_master();
master_ptr0 = master_ptr;
if (master_ptr0 == master_ptr->s->primary)
@ -2146,6 +2148,7 @@ mb_for_species_aq(int n)
* by coef, usually moles.
* mb_unknowns.coef - coefficient of s[n] in equation or relation
*/
//int i, j;
int i;
class master *master_ptr;
class unknown *unknown_ptr;
@ -2222,6 +2225,7 @@ mb_for_species_aq(int n)
*/
if (use.Get_surface_ptr() != NULL && s[n]->type < H2O && dl_type_x != cxxSurface::NO_DL)
{
//j = 0;
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type == SURFACE_CB)
@ -2233,6 +2237,7 @@ mb_for_species_aq(int n)
store_mb_unknowns(unknown_ptr, s_diff_layer[n][charge_ptr->Get_name()].Get_g_moles_address(),
s[n]->z, s_diff_layer[n][charge_ptr->Get_name()].Get_dg_g_moles_address());
//j++;
}
}
}

12
print.cpp
View File

@ -339,8 +339,12 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
{
LDBLE exp_g = charge_ptr->Get_g_map()[1].Get_g() * mass_water_aq_x / mass_water_surface + 1;
LDBLE psi_DL = -log(exp_g) * R_KJ_DEG_MOL * tk_x / F_KJ_V_EQ;
output_msg(sformatf(
"\n\tTotal moles in diffuse layer (excluding water), Donnan calculation."));
if (use.Get_surface_ptr()->Get_correct_GC())
output_msg(sformatf(
"\n\tTotal moles in diffuse layer (excluding water), Donnan corrected to match Poisson-Boltzmann."));
else
output_msg(sformatf(
"\n\tTotal moles in diffuse layer (excluding water), Donnan calculation."));
output_msg(sformatf(
"\n\tDonnan Layer potential, psi_DL = %10.3e V.\n\tBoltzmann factor, exp(-psi_DL * F / RT) = %9.3e (= c_DL / c_free if z is +1).\n\n",
psi_DL, exp_g));
@ -2234,7 +2238,7 @@ print_totals(void)
(double) calc_solution_volume()));
}
/* VP: Density End */
#ifdef NPP
//#ifdef NPP
if (print_viscosity)
{
output_msg(sformatf("%45s%9.5f", "Viscosity (mPa s) = ",
@ -2251,7 +2255,7 @@ print_totals(void)
}
else output_msg(sformatf("\n"));
}
#endif
//#endif
output_msg(sformatf("%45s%7.3f\n", "Activity of water = ",
exp(s_h2o->la * LOG_10)));
output_msg(sformatf("%45s%11.3e\n", "Ionic strength (mol/kgw) = ",

51
read.cpp
View File

@ -6283,7 +6283,7 @@ read_surface(void)
* ERROR if error occurred reading data
*
*/
int n_user;
int n_user, i1;
LDBLE conc;
const char* cptr, *cptr1;
std::string token, token1, name;
@ -6306,9 +6306,10 @@ read_surface(void)
"ccm", /* 13 */
"equilibrium", /* 14 */
"site_units", /* 15 */
"ddl" /* 16 */
"ddl", /* 16 */
"donnan_factors" /* 17 */
};
int count_opt_list = 17;
int count_opt_list = 18;
/*
* kin_surf is for Surfaces, related to kinetically reacting minerals
* they are defined if "sites" is followed by mineral name:
@ -6413,7 +6414,7 @@ read_surface(void)
if (thickness != 0)
{
error_msg
("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8",
("You must enter EITHER thickness OR Debye lengths (1/k),\n and relative DDL viscosity, DDL limit.\nCorrect is (for example): -donnan 1e-8 viscosity 0.5 limit 0.9 correct_GC true\n or (default values): -donnan debye_lengths 1 viscosity 1 limit 0.8 correct_GC false",
CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
@ -6436,6 +6437,23 @@ read_surface(void)
break;
}
}
else if (token[0] == 'C' || token[0] == 'c')
{
copy_token(token1, &next_char);
if (token1[0] == 'T' || token1[0] == 't' || token1[0] == 'F' || token1[0] == 'f')
{
temp_surface.Set_correct_GC(get_true_false(token1.c_str(), TRUE) == TRUE);
continue;
} else
{
error_msg
("Expected True or False for correct_GC (which brings co-ion concentrations closer to their integrated double layer value).",
CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
break;
}
}
else if (token[0] == 'V' || token[0] == 'v')
{
int j = copy_token(token1, &next_char);
@ -6570,6 +6588,31 @@ read_surface(void)
case 16: /* ddl */
temp_surface.Set_type(cxxSurface::DDL);
break;
case 17: /* Donnan_factors */
temp_surface.Donnan_factors.clear();
i1 = 0;
for (;;)
{
int i = copy_token(token, &next_char);
if (i == DIGIT && i1 < 8)
{
(void)sscanf(token.c_str(), SCANFORMAT, &dummy);
temp_surface.Donnan_factors.push_back(dummy);
i1++;
continue;
}
else if (i != EMPTY || i1 > 8)
{
error_msg
("Expected at most 8 numbers for the Donnan_factors for co- and counter-ions,\n z *= cgc[0] * (mu_x**(cgc[1] * nDbl**cgc[2] * (abs(surf_chrg_eq / A_surf / 1e-6)**cgc[3] * mu_x**(cgc[4])",
CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
break;
}
break;
}
break;
case OPTION_DEFAULT:
/*
* Read surface component

3
spread.cpp
View File

@ -1141,6 +1141,7 @@ copy_token_tab(std::string& token, const char **cptr)
* EOL,
* UNKNOWN.
*/
//int i, return_value;
int return_value;
char c;
/*
@ -1180,6 +1181,7 @@ copy_token_tab(std::string& token, const char **cptr)
/*
* Begin copying to token
*/
//i = 0;
for (;;)
{
c = **cptr;
@ -1196,6 +1198,7 @@ copy_token_tab(std::string& token, const char **cptr)
{
token.push_back(c);
(*cptr)++;
//i++;
}
}
return (return_value);

3
sundialsmath.cpp
View File

@ -59,7 +59,8 @@
#include <stdio.h>
#include <math.h>
//#include <math.h>
#include <cmath>
#include "sundialsmath.h"
#include "sundialstypes.h"

7
tally.cpp
View File

@ -798,6 +798,7 @@ build_tally_table(void)
*/
int j, k, l, p, save_print_use;
size_t n;
//int count_tt_pure_phase, count_tt_ss_phase, count_tt_kinetics;
class phase *phase_ptr;
char token[MAX_LENGTH];
const char* cptr;
@ -870,6 +871,7 @@ build_tally_table(void)
/*
* Count pure phases
*/
//count_tt_pure_phase = 0;
if (Rxn_pp_assemblage_map.size() > 0)
{
/*
@ -902,6 +904,7 @@ build_tally_table(void)
/*
* Add to table
*/
//count_tt_pure_phase++;
n = count_tally_table_columns;
extend_tally_table();
tally_table[n].name = phase_ptr->name;
@ -928,6 +931,7 @@ build_tally_table(void)
/*
* Add solid-solution pure phases
*/
//count_tt_ss_phase = 0;
if (Rxn_ss_assemblage_map.size() > 0)
{
/*
@ -960,6 +964,7 @@ build_tally_table(void)
/*
* Add to table
*/
//count_tt_ss_phase++;
n = count_tally_table_columns;
extend_tally_table();
tally_table[n].name = phase_ptr->name;
@ -977,6 +982,7 @@ build_tally_table(void)
/*
* Add kinetic reactants
*/
//count_tt_kinetics = 0;
if (Rxn_kinetics_map.size() > 0)
{
std::map<int, cxxKinetics>::iterator it;
@ -1000,6 +1006,7 @@ build_tally_table(void)
/*
* Add to table
*/
//count_tt_kinetics++;
n = count_tally_table_columns;
extend_tally_table();
tally_table[n].name = string_hsave(kinetics_comp_ptr->Get_rate_name().c_str());

12
tidy.cpp
View File

@ -982,7 +982,8 @@ tidy_gas_phase(void)
#ifdef NPP
if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#else
if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
//if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#endif
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
@ -1015,7 +1016,8 @@ tidy_gas_phase(void)
#ifdef NPP
if (!isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#else
if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
//if (gas_phase_ptr->Get_gas_comps()[j].Get_p_read() != NAN)
if (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#endif
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
@ -1690,7 +1692,8 @@ tidy_ss_assemblage(void)
#ifdef NPP
if (isnan(comp_ptr->Get_moles()))
#else
if (comp_ptr->Get_moles() == NAN)
//if (comp_ptr->Get_moles() == NAN)
if (std::isnan(comp_ptr->Get_moles()))
#endif
{
input_error++;
@ -3021,7 +3024,8 @@ tidy_isotopes(void)
#ifdef NPP
if (!isnan(master[k]->isotope_ratio_uncertainty))
#else
if (master[k]->isotope_ratio_uncertainty != NAN)
//if (master[k]->isotope_ratio_uncertainty != NAN)
if (!std::isnan(master[k]->isotope_ratio_uncertainty))
#endif
{
temp_iso.Set_ratio_uncertainty_defined(true);

54
transport.cpp
View File

@ -1888,16 +1888,17 @@ fill_spec(int l_cell_no, int ref_cell)
if (por_il < interlayer_Dpor_lim)
por_il = viscos_il_f = 0.0;
/*
* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos
* modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959
*/
viscos = viscos_0;
/*
* put temperature factor in por_factor which corrects for porous medium...
*/
viscos_f *= tk_x * viscos_0_25 / (298.15 * viscos);
viscos_il_f *= tk_x * viscos_0_25 / (298.15 * viscos);
sol_D[l_cell_no].viscos_f = tk_x * viscos_0_25 / (298.15 * viscos);
dum = viscos_0_25 / viscos;
viscos_f *= dum;
viscos_il_f *= dum;
sol_D[l_cell_no].viscos_f = dum;
count_spec = count_exch_spec = 0;
/*
@ -1908,6 +1909,8 @@ fill_spec(int l_cell_no, int ref_cell)
qsort(&species_list[0], species_list.size(),
sizeof(class species_list), sort_species_name);
}
if (correct_Dw)
calc_SC();
for (i = 0; i < (int)species_list.size(); i++)
{
@ -2097,7 +2100,7 @@ fill_spec(int l_cell_no, int ref_cell)
}
if (correct_Dw)
{
calc_SC(); // note that neutral species are corrected as if z = 1, but is viscosity-dependent
//calc_SC(); // removed that neutral species are corrected as if z = 1, but is viscosity-dependent
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr * viscos_f;
}
if (l_cell_no <= count_cells + 1 && sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn) > diffc_max)
@ -2962,20 +2965,20 @@ diffuse_implicit(LDBLE DDt, int stagnant)
if (!strcmp(ct[icell].m_s[cp].name, "H"))
{
dummy = ct[icell].m_s[cp].tot1;
if (dV_dcell || (icell > 0 && icell <= ilast))
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
sptr2->Set_total_h(sptr2->Get_total_h() + dummy);
if (sptr_stag)
{
dummy = ct[icell].m_s[cp].tot_stag;
if (dV_dcell || (icell > 0 && icell <= ilast))
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
if (icell == c)
{
// mix the boundary solution with the stagnant column end-cell
dummy += ct[icell + 1].m_s[cp].tot_stag;
if (dV_dcell)
if (dV_dcell || fix_current)
sptr2->Set_total_h(sptr2->Get_total_h() - ct[icell + 1].m_s[cp].tot_stag);
}
sptr_stag->Set_total_h(sptr_stag->Get_total_h() + dummy);
@ -2985,19 +2988,19 @@ diffuse_implicit(LDBLE DDt, int stagnant)
if (!strcmp(ct[icell].m_s[cp].name, "O"))
{
dummy = ct[icell].m_s[cp].tot1;
if (dV_dcell || (icell > 0 && icell <= ilast))
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
sptr2->Set_total_o(sptr2->Get_total_o() + dummy);
if (sptr_stag)
{
dummy = ct[icell].m_s[cp].tot_stag;
if (dV_dcell || (icell > 0 && icell <= ilast))
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
if (icell == c)
{
dummy += ct[icell + 1].m_s[cp].tot_stag;
if (dV_dcell)
if (dV_dcell || fix_current)
sptr2->Set_total_o(sptr2->Get_total_o() - ct[icell + 1].m_s[cp].tot_stag);
}
sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy);
@ -3018,7 +3021,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
// 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)))
(dV_dcell || fix_current || (icell > 0 && icell <= ilast)))
{
dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name);
if (ct[icell].dl_s > 1e-8)
@ -3039,7 +3042,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 || fix_current || (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)
@ -3057,7 +3060,7 @@ diffuse_implicit(LDBLE DDt, int stagnant)
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 (dV_dcell || (icell > 0 && icell <= ilast))
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
{
dum = ct[icell].m_s[cp].tot1;
if (stagnant)
@ -3078,10 +3081,10 @@ 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 || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
{
dum = ct[icell].m_s[cp].tot1;
if (stagnant && icell == c && dV_dcell)
if (stagnant && icell == c && (dV_dcell || fix_current))
dum -= ct[c1].m_s[cp].tot_stag;
dum2 += dum;
sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol);
@ -3134,7 +3137,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 (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)
if ((dV_dcell/* || fix_current*/) && 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)
@ -4828,10 +4831,10 @@ Step (from cell 1 to count_cells + 1):
cxxSurface *surface_ptr1, *surface_ptr2;
LDBLE viscos_f;
/*
* temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
* temperature and viscosity correction for MCD coefficient, D_T = D_298 * viscos_298 / viscos
*/
viscos_f = viscos_0;
viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f);
viscos_f = viscos_0_25 / viscos_f;
//n1 = 0;
//n2 = n1 + 1;
@ -5502,7 +5505,7 @@ diff_stag_surf(int mobile_cell)
* temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
*/
viscos_f = viscos_0;
viscos_f = tk_x * viscos_0_25 / (298.15 * viscos_f);
viscos_f = viscos_0_25 / viscos_f;
cxxSurface surface_n1, surface_n2;
cxxSurface *surface_n1_ptr = &surface_n1;
@ -5991,8 +5994,7 @@ viscosity(void)
}
// find B * m and D * m * mu^d3
Bc += (s_x[i]->Jones_Dole[0] +
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) *
t1;
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
// define f_I from the exponent of the D * m^d3 term...
if (s_x[i]->Jones_Dole[5] >= 1)
t2 = mu_x / 3 / s_x[i]->Jones_Dole[5];
@ -6074,8 +6076,8 @@ viscosity(void)
viscos += (viscos_0 * (Bc + Dc));
if (viscos < 0)
{
viscos = 0;
warning_msg("viscosity < 0, reset to 0.");
viscos = viscos_0;
warning_msg("viscosity < 0, reset to viscosity of pure water");
}
return viscos;