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 '311fa04fe162e939e01419a658681ec819a534ea'

Browse Source
This commit is contained in:
Scott R Charlton 2018-08-01 18:27:29 -06:00
commit 80462a2d42
22 changed files with 516 additions and 354 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/phreeqcpp/PBasic.cpp
View File

@ -1655,6 +1655,9 @@ listtokens(FILE * f, tokenrec * l_buf)
case tokt_sc:
output_msg("T_SC");
break;
case tokiterations:
output_msg("ITERATIONS");
break;
}
l_buf = l_buf->next;
}
@ -3710,6 +3713,9 @@ factor(struct LOC_exec * LINK)
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_t_sc(str);
}
break;
case tokiterations:
n.UU.val = (parse_all) ? 0 : PhreeqcPtr->overall_iterations;
break;
case toklog10:
{
@ -7333,9 +7339,10 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("phase_vm", PBasic::tokphase_vm),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("current_a", PBasic::tokcurrent_a),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("pot_v", PBasic::tokpot_v),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("t_sc", PBasic::tokt_sc),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("t_sc", PBasic::tokt_sc),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("setdiff_c", PBasic::toksetdiff_c),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("aphi", PBasic::tokaphi)
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("aphi", PBasic::tokaphi),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("iterations", PBasic::tokiterations)
};
std::map<const std::string, PBasic::BASIC_TOKEN> PBasic::command_tokens(temp_tokens, temp_tokens + sizeof temp_tokens / sizeof temp_tokens[0]);

3
src/phreeqcpp/PBasic.h
View File

@ -336,7 +336,8 @@ public:
tokcurrent_a,
tokpot_v,
tokt_sc,
tokaphi
tokaphi,
tokiterations
};
#if !defined(PHREEQCI_GUI)

7
src/phreeqcpp/Phreeqc.cpp
View File

@ -695,7 +695,7 @@ void Phreeqc::init(void)
stop_program = FALSE;
incremental_reactions = FALSE;
count_strings = 0;
array = NULL;
my_array = NULL;
delta = NULL;
residual = NULL;
input_error = 0;
@ -705,6 +705,7 @@ void Phreeqc::init(void)
iterations = 0;
gamma_iterations = 0;
run_reactions_iterations= 0;
overall_iterations = 0;
max_line = MAX_LINE;
line = NULL;
line_save = NULL;
@ -712,6 +713,8 @@ void Phreeqc::init(void)
debug_model = FALSE;
debug_prep = FALSE;
debug_set = FALSE;
debug_mass_action = FALSE;
debug_mass_balance = FALSE;
debug_diffuse_layer = FALSE;
debug_inverse = FALSE;
#ifdef USE_LONG_DOUBLE
@ -1921,7 +1924,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
max_line = MAX_LINE;
line = NULL;
line_save = NULL;
LOG_10 = log(10.0);
LOG_10 = LOG_10;
debug_model = FALSE;
debug_prep = FALSE;
debug_set = FALSE;

5
src/phreeqcpp/Phreeqc.h
View File

@ -1625,7 +1625,7 @@ protected:
int count_strings;
int max_strings;
LDBLE *array;
LDBLE *my_array;
LDBLE *delta;
LDBLE *residual;
@ -1637,6 +1637,7 @@ protected:
int iterations;
int gamma_iterations;
int run_reactions_iterations;
int overall_iterations;
int max_line;
char *line;
@ -1646,6 +1647,8 @@ protected:
int debug_model;
int debug_prep;
int debug_mass_action;
int debug_mass_balance;
int debug_set;
int debug_diffuse_layer;
int debug_inverse;

2
src/phreeqcpp/SelectedOutput.cpp
View File

@ -38,6 +38,7 @@ SelectedOutput::SelectedOutput(int n, PHRQ_io *io)
this->water = false;
this->charge_balance = false;
this->percent_error = false;
this->new_line = true;
// as-is set flags
//
@ -62,6 +63,7 @@ SelectedOutput::SelectedOutput(int n, PHRQ_io *io)
this->set_water = false;
this->set_charge_balance = false;
this->set_percent_error = false;
this->set_new_line = false;
}
SelectedOutput::~SelectedOutput(void)

7
src/phreeqcpp/SelectedOutput.h
View File

@ -82,6 +82,7 @@ public:
inline bool Get_water(void)const {return this->water;}
inline bool Get_charge_balance(void)const {return this->charge_balance;}
inline bool Get_percent_error(void)const {return this->percent_error;}
inline bool Get_new_line(void)const {return this->new_line; }
// as-is setters
inline void Set_user_punch(bool tf) {this->user_punch = tf; this->set_user_punch = true;}
@ -105,6 +106,7 @@ public:
inline void Set_water(bool tf) {this->water = tf; this->set_water = true;}
inline void Set_charge_balance(bool tf) {this->charge_balance = tf; this->set_charge_balance = true;}
inline void Set_percent_error(bool tf) {this->percent_error = tf; this->set_percent_error = true;}
inline void Set_new_line(bool tf) {this->new_line = tf; this->set_new_line = true;}
// set flag getters
inline bool was_set_user_punch()const {return this->set_user_punch;}
@ -128,7 +130,8 @@ public:
inline bool was_set_water()const {return this->set_water;}
inline bool was_set_charge_balance()const {return this->set_charge_balance;}
inline bool was_set_percent_error()const {return this->set_percent_error;}
inline bool was_set_new_line()const {return this->set_new_line;}
protected:
// vectors
@ -177,6 +180,7 @@ protected:
bool water;
bool charge_balance;
bool percent_error;
bool new_line;
// as-is set flags
bool set_user_punch;
@ -200,5 +204,6 @@ protected:
bool set_water;
bool set_charge_balance;
bool set_percent_error;
bool set_new_line;
};
#endif // !defined(SELECTEDOUTPUT_H_INCLUDED)

10
src/phreeqcpp/class_main.cpp
View File

@ -277,8 +277,8 @@ write_banner(void)
/* version */
#ifdef NPP
//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.4 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.4.4");
//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.6 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.4.6");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -302,7 +302,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "April 5, 2018");
len = sprintf(buffer, "%s", "June 20, 2018");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -492,12 +492,12 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
screen_msg(sformatf("Database file: %s\n\n", token));
strcpy(db_file, token);
#ifdef NPP
//output_msg(sformatf("Using PHREEQC: version 3.4.4 Ampère, compiled April 5, 2018\n"));
//output_msg(sformatf("Using PHREEQC: version 3.4.6 Ampère, compiled June 20, 2018\n"));
#endif
output_msg(sformatf(" Input file: %s\n", in_file));
output_msg(sformatf(" Output file: %s\n", out_file));
#ifdef NPP
output_msg(sformatf("Using PHREEQC: version 3.4.4, compiled April 5, 2018\n"));
output_msg(sformatf("Using PHREEQC: version 3.4.6, compiled June 20, 2018\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token));
#ifdef NPP

18
src/phreeqcpp/gases.cpp
View File

@ -87,7 +87,7 @@ build_fixed_volume_gas(void)
*/
if (debug_prep == TRUE)
{
output_msg(sformatf( "\n\tMass balance summations %s.\n\n",
output_msg(sformatf( "\n\tMass balance summations %s.\n",
phase_ptr->name));
}
@ -211,27 +211,27 @@ build_fixed_volume_gas(void)
}
col = master_ptr->unknown->number;
coef = coef_elt * rxn_ptr->coef;
store_jacob(&(gas_unknowns[i]->moles),
&(array[row + col]), coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
master_ptr->s->name, (double) coef,
row / (count_unknowns + 1), col));
}
store_jacob(&(gas_unknowns[i]->moles),
&(my_array[row + col]), coef);
}
if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE)
{
/* derivative wrt total moles of gas */
store_jacob(&(phase_ptr->fraction_x),
&(array[row + gas_unknown->number]), coef_elt);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
"gas moles", (double) elt_list[j].coef,
row / (count_unknowns + 1),
gas_unknown->number));
}
store_jacob(&(phase_ptr->fraction_x),
&(my_array[row + gas_unknown->number]), coef_elt);
}
}
/*
@ -303,13 +303,13 @@ build_fixed_volume_gas(void)
}
col = master_ptr->unknown->number;
coef = rxn_ptr->coef;
store_jacob(&(phase_ptr->p_soln_x), &(array[row + col]), coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
master_ptr->s->name, (double) coef,
row / (count_unknowns + 1), col));
}
store_jacob(&(phase_ptr->p_soln_x), &(my_array[row + col]), coef);
}
}
}

18
src/phreeqcpp/integrate.cpp
View File

@ -604,8 +604,7 @@ initial_surface_water(void)
if (rd > rd_limit)
{
mass_water_surfaces_x =
use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 -
ddl_limit);
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;
rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
@ -614,15 +613,13 @@ initial_surface_water(void)
}
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)
continue;
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
l_s =charge_ptr->Get_specific_area() *
charge_ptr->Get_grams();
l_s = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
charge_ptr->Set_mass_water(mass_water_surfaces_x * l_s / sum_surfs);
}
}
@ -754,7 +751,8 @@ calc_all_donnan(void)
/*
* calculate g for each surface...
*/
initial_surface_water();
if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths()) // DL_pitz
initial_surface_water();
converge = TRUE;
for (int j = 0; j < count_unknowns; j++)
{
@ -835,9 +833,9 @@ calc_all_donnan(void)
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 < 0 && z < 0)
|| (surf_chrg_eq > 0 && z > 0)))
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;
if (new_g2 <= -ratio_aq)
new_g2 = -ratio_aq + G_TOL * 1e-3;

170
src/phreeqcpp/inverse.cpp
View File

@ -266,11 +266,11 @@ setup_inverse(struct inverse *inv_ptr)
/*
* Malloc space for arrays
*/
array = (LDBLE *) free_check_null(array);
array =
my_array = (LDBLE *) free_check_null(my_array);
my_array =
(LDBLE *) PHRQ_malloc((size_t) max_column_count * max_row_count *
sizeof(LDBLE));
if (array == NULL)
if (my_array == NULL)
malloc_error();
array1 =
@ -349,7 +349,7 @@ setup_inverse(struct inverse *inv_ptr)
(size_t) max_column_count * sizeof(LDBLE));
for (i = 0; i < max_row_count; i++)
{
memcpy((void *) &(array[i * max_column_count]), (void *) &(inv_zero[0]),
memcpy((void *) &(my_array[i * max_column_count]), (void *) &(inv_zero[0]),
(size_t) max_column_count * sizeof(LDBLE));
}
/*
@ -443,11 +443,11 @@ setup_inverse(struct inverse *inv_ptr)
{
if (master[j]->in >= 0)
{
array[master[j]->in * max_column_count + i] =
my_array[master[j]->in * max_column_count + i] =
f * master[j]->total;
if (master[j]->s == s_eminus)
{
array[master[j]->in * max_column_count + i] = 0.0;
my_array[master[j]->in * max_column_count + i] = 0.0;
}
}
}
@ -474,7 +474,7 @@ setup_inverse(struct inverse *inv_ptr)
}
if (fabs(cb) < toler)
cb = 0.0;
array[(row_charge + i) * max_column_count + i] = cb;
my_array[(row_charge + i) * max_column_count + i] = cb;
}
/* mass_balance: phase data */
@ -520,11 +520,11 @@ setup_inverse(struct inverse *inv_ptr)
coef = master_ptr->coef;
if (coef <= 0)
coef = 1.0;
array[row * max_column_count + column] =
my_array[row * max_column_count + column] =
rxn_ptr->token[j].coef * coef;
}
row = master_alk->in; /* include alkalinity for phase */
array[row * max_column_count + column] = calc_alk(rxn_ptr);
my_array[row * max_column_count + column] = calc_alk(rxn_ptr);
}
/* mass balance: redox reaction data */
@ -572,14 +572,14 @@ setup_inverse(struct inverse *inv_ptr)
assert(row * max_column_count + column < max_column_count * max_row_count);
assert(row >= 0);
assert(column >= 0);
array[row * max_column_count + column] =
my_array[row * max_column_count + column] =
rxn_ptr->token[j].coef;
/* if coefficient of element is not 1.0 in master species */
if (j != 0)
array[row * max_column_count + column] /= coef;
my_array[row * max_column_count + column] /= coef;
}
row = master_alk->in; /* include alkalinity for redox reaction */
array[row * max_column_count + column] =
my_array[row * max_column_count + column] =
(calc_alk(rxn_ptr) - inv_ptr->elts[i].master->s->alk) / coef;
}
}
@ -594,15 +594,15 @@ setup_inverse(struct inverse *inv_ptr)
{
if (j < (inv_ptr->count_solns - 1))
{
array[row * max_column_count + column] = 1.0;
my_array[row * max_column_count + column] = 1.0;
}
else
{
array[row * max_column_count + column] = -1.0;
my_array[row * max_column_count + column] = -1.0;
}
if (inv_ptr->elts[i].master->s == s_eminus)
{
array[row * max_column_count + column] = 0.0;
my_array[row * max_column_count + column] = 0.0;
}
sprintf(token, "%s %d", row_name[row], j);
col_name[column] = string_hsave(token);
@ -664,19 +664,19 @@ setup_inverse(struct inverse *inv_ptr)
solution_ptr = Utilities::Rxn_find(Rxn_solution_map, inv_ptr->solns[i]);
if (i < inv_ptr->count_solns - 1)
{
array[count_rows * max_column_count + i] =
my_array[count_rows * max_column_count + i] =
1.0 / gfw_water * solution_ptr->Get_mass_water();
}
else
{
array[count_rows * max_column_count + inv_ptr->count_solns - 1] =
my_array[count_rows * max_column_count + inv_ptr->count_solns - 1] =
-1.0 / gfw_water * solution_ptr->Get_mass_water();
}
}
/* coefficient for water uncertainty */
if (inv_ptr->water_uncertainty > 0)
{
array[count_rows * max_column_count + col_water] = 1.0;
my_array[count_rows * max_column_count + col_water] = 1.0;
}
row_name[count_rows] = string_hsave("H2O");
row_water = count_rows;
@ -686,8 +686,8 @@ setup_inverse(struct inverse *inv_ptr)
* Final solution fraction equals 1.0
*/
array[count_rows * max_column_count + inv_ptr->count_solns - 1] = 1.0;
array[count_rows * max_column_count + count_unknowns] = 1.0;
my_array[count_rows * max_column_count + inv_ptr->count_solns - 1] = 1.0;
my_array[count_rows * max_column_count + count_unknowns] = 1.0;
row_name[count_rows] = string_hsave("fract, final");
count_rows++;
@ -709,10 +709,10 @@ setup_inverse(struct inverse *inv_ptr)
{
coef = -1.0;
}
array[count_rows * max_column_count + column] = coef;
my_array[count_rows * max_column_count + column] = coef;
if (inv_ptr->elts[j].master->s == s_eminus)
{
array[count_rows * max_column_count + column] = 0.0;
my_array[count_rows * max_column_count + column] = 0.0;
}
}
sprintf(token, "%s %d", "charge", i);
@ -730,12 +730,12 @@ setup_inverse(struct inverse *inv_ptr)
if (inv_ptr->dalk_dph[i] != 0 || inv_ptr->dalk_dc[i] != 0)
{
column = col_ph + i;
array[count_rows * max_column_count + column] =
my_array[count_rows * max_column_count + column] =
inv_ptr->dalk_dph[i];
column = col_epsilon + i_alk * inv_ptr->count_solns + i;
array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + column] = -1.0;
column = col_epsilon + i_carb * inv_ptr->count_solns + i;
array[count_rows * max_column_count + column] =
my_array[count_rows * max_column_count + column] =
inv_ptr->dalk_dc[i];
}
sprintf(token, "%s %d", "dAlk", i);
@ -782,7 +782,7 @@ setup_inverse(struct inverse *inv_ptr)
else
{
coef =
array[inv_ptr->elts[j].master->in * max_column_count +
my_array[inv_ptr->elts[j].master->in * max_column_count +
i] * coef;
coef = fabs(coef);
}
@ -795,7 +795,7 @@ setup_inverse(struct inverse *inv_ptr)
{
for (k = 0; k < count_rows; k++)
{
array[k * max_column_count + column] = 0.0;
my_array[k * max_column_count + column] = 0.0;
}
continue;
}
@ -807,12 +807,12 @@ setup_inverse(struct inverse *inv_ptr)
if (coef < toler)
{
array[(column - col_epsilon) * max_column_count + column] =
my_array[(column - col_epsilon) * max_column_count + column] =
SCALE_EPSILON / toler;
}
else
{
array[(column - col_epsilon) * max_column_count + column] =
my_array[(column - col_epsilon) * max_column_count + column] =
SCALE_EPSILON / coef;
}
@ -826,15 +826,15 @@ setup_inverse(struct inverse *inv_ptr)
{
f = 1.0;
}
array[count_rows * max_column_count + column] = 1.0 * f;
array[count_rows * max_column_count + i] = -coef * f;
my_array[count_rows * max_column_count + column] = 1.0 * f;
my_array[count_rows * max_column_count + i] = -coef * f;
sprintf(token, "%s %s", inv_ptr->elts[j].master->elt->name,
"eps+");
row_name[count_rows] = string_hsave(token);
count_rows++;
/* set lower limit of change in negative direction */
conc = array[inv_ptr->elts[j].master->in * max_column_count + i];
conc = my_array[inv_ptr->elts[j].master->in * max_column_count + i];
/* if concentration is zero, only positive direction allowed */
if (conc == 0.0)
@ -861,8 +861,8 @@ setup_inverse(struct inverse *inv_ptr)
inv_ptr->elts[j].master->elt->name))
coef = fabs(conc) + toler;
array[count_rows * max_column_count + i] = -coef * f;
array[count_rows * max_column_count + column] = -1.0 * f;
my_array[count_rows * max_column_count + i] = -coef * f;
my_array[count_rows * max_column_count + column] = -1.0 * f;
sprintf(token, "%s %s", inv_ptr->elts[j].master->elt->name,
"eps-");
row_name[count_rows] = string_hsave(token);
@ -882,21 +882,21 @@ setup_inverse(struct inverse *inv_ptr)
/* scale epsilon in optimization equation */
array[(column - col_epsilon) * max_column_count + column] =
my_array[(column - col_epsilon) * max_column_count + column] =
SCALE_EPSILON / coef;
/* set upper limit of change in positive direction */
array[count_rows * max_column_count + column] = 1.0;
array[count_rows * max_column_count + i] = -coef;
my_array[count_rows * max_column_count + column] = 1.0;
my_array[count_rows * max_column_count + i] = -coef;
sprintf(token, "%s %s", "pH", "eps+");
row_name[count_rows] = string_hsave(token);
count_rows++;
/* set lower limit of change in negative direction */
array[count_rows * max_column_count + column] = -1.0;
array[count_rows * max_column_count + i] = -coef;
my_array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + i] = -coef;
sprintf(token, "%s %s", "pH", "eps-");
row_name[count_rows] = string_hsave(token);
count_rows++;
@ -911,16 +911,16 @@ setup_inverse(struct inverse *inv_ptr)
if (coef > 0.0)
{
/* set upper limit of change in positive direction */
array[count_rows * max_column_count + column] = 1.0;
array[count_rows * max_column_count + count_unknowns] = coef;
my_array[count_rows * max_column_count + column] = 1.0;
my_array[count_rows * max_column_count + count_unknowns] = coef;
sprintf(token, "%s %s", "water", "eps+");
row_name[count_rows] = string_hsave(token);
count_rows++;
/* set lower limit of change in negative direction */
array[count_rows * max_column_count + column] = -1.0;
array[count_rows * max_column_count + count_unknowns] = coef;
my_array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + count_unknowns] = coef;
sprintf(token, "%s %s", "water", "eps-");
row_name[count_rows] = string_hsave(token);
count_rows++;
@ -952,12 +952,12 @@ setup_inverse(struct inverse *inv_ptr)
/* scale epsilon in optimization equation */
array[(column - col_epsilon) * max_column_count +
my_array[(column - col_epsilon) * max_column_count +
column] = SCALE_EPSILON / coef;
/* set upper limit of change in positive direction */
array[count_rows * max_column_count + column] = 1.0;
array[count_rows * max_column_count + i] = -coef;
my_array[count_rows * max_column_count + column] = 1.0;
my_array[count_rows * max_column_count + i] = -coef;
sprintf(token, "%d%s %s",
(int) kit->second.Get_isotope_number(),
kit->second.Get_elt_name().c_str(), "eps+");
@ -966,8 +966,8 @@ setup_inverse(struct inverse *inv_ptr)
/* set lower limit of change in negative direction */
array[count_rows * max_column_count + column] = -1.0;
array[count_rows * max_column_count + i] = -coef;
my_array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + i] = -coef;
sprintf(token, "%d%s %s",
(int) kit->second.Get_isotope_number(),
kit->second.Get_elt_name().c_str(), "eps-");
@ -1012,7 +1012,7 @@ setup_inverse(struct inverse *inv_ptr)
*/
for (i = 0; i < max_column_count; i++)
{
array[row_water * max_column_count + i] *= SCALE_WATER;
my_array[row_water * max_column_count + i] *= SCALE_WATER;
}
/*
* Arrays are complete
@ -1035,7 +1035,7 @@ setup_inverse(struct inverse *inv_ptr)
k = 0;
}
output_msg(sformatf( "%11.2e",
(double) array[i * max_column_count + j]));
(double) my_array[i * max_column_count + j]));
k++;
}
if (k != 0)
@ -1332,7 +1332,7 @@ solve_inverse(struct inverse *inv_ptr)
output_msg(sformatf( "\tNumber of calls to cl1: %d\n",
count_calls));
}
array = (LDBLE *) free_check_null(array);
my_array = (LDBLE *) free_check_null(my_array);
delta = (LDBLE *) free_check_null(delta);
array1 = (LDBLE *) free_check_null(array1);
inv_zero = (LDBLE *) free_check_null(inv_zero);
@ -1462,7 +1462,7 @@ solve_with_mask(struct inverse *inv_ptr, unsigned long cur_bits)
memcpy((void *) &(delta_save[0]), (void *) &(inv_zero[0]),
(size_t) max_column_count * sizeof(LDBLE));
shrink(inv_ptr, array, array1,
shrink(inv_ptr, my_array, array1,
&k, &l, &m, &n, cur_bits, delta2, col_back, row_back);
/*
* Save delta constraints
@ -2523,7 +2523,7 @@ range(struct inverse *inv_ptr, unsigned long cur_bits)
/*
* Copy equations
*/
memcpy((void *) &(array1[0]), (void *) &(array[0]),
memcpy((void *) &(array1[0]), (void *) &(my_array[0]),
(size_t) max_column_count * max_row_count * sizeof(LDBLE));
memcpy((void *) &(delta2[0]), (void *) &(delta[0]),
(size_t) max_column_count * sizeof(LDBLE));
@ -2971,7 +2971,7 @@ check_solns(struct inverse *inv_ptr)
/*
* Copy equations
*/
memcpy((void *) &(array1[0]), (void *) &(array[0]),
memcpy((void *) &(array1[0]), (void *) &(my_array[0]),
(size_t) max_column_count * max_row_count * sizeof(LDBLE));
memcpy((void *) &(delta2[0]), (void *) &(delta[0]),
(size_t) max_column_count * sizeof(LDBLE));
@ -3114,16 +3114,16 @@ post_mortem(void)
sum = 0;
for (j = 0; j < count_unknowns; j++)
{
sum += inv_delta1[j] * array[i * max_column_count + j];
sum += inv_delta1[j] * my_array[i * max_column_count + j];
}
if (equal(sum, array[(i * max_column_count) + count_unknowns], toler)
if (equal(sum, my_array[(i * max_column_count) + count_unknowns], toler)
== FALSE)
{
output_msg(sformatf(
"\tERROR: equality not satisfied for %s, %e.\n",
row_name[i],
(double) (sum - array[(i * max_column_count) + count_unknowns])));
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
}
}
/*
@ -3134,15 +3134,15 @@ post_mortem(void)
sum = 0;
for (j = 0; j < count_unknowns; j++)
{
sum += inv_delta1[j] * array[i * max_column_count + j];
sum += inv_delta1[j] * my_array[i * max_column_count + j];
}
if (sum > array[(i * max_column_count) + count_unknowns] + toler)
if (sum > my_array[(i * max_column_count) + count_unknowns] + toler)
{
output_msg(sformatf(
"\tERROR: inequality not satisfied for %s, %e\n",
row_name[i],
(double) (sum - array[(i * max_column_count) + count_unknowns])));
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
}
}
/*
@ -3191,15 +3191,15 @@ test_cl1_solution(void)
sum = 0;
for (j = 0; j < count_unknowns; j++)
{
sum += inv_delta1[j] * array[i * max_column_count + j];
sum += inv_delta1[j] * my_array[i * max_column_count + j];
}
if (equal(sum, array[(i * max_column_count) + count_unknowns], toler) == FALSE)
if (equal(sum, my_array[(i * max_column_count) + count_unknowns], toler) == FALSE)
{
if (debug_inverse)
{
output_msg(sformatf("\tERROR: equality not satisfied for %s, %e.\n", row_name[i],
(double) (sum - array[(i * max_column_count) + count_unknowns])));
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
}
rv = false;
}
@ -3212,17 +3212,17 @@ test_cl1_solution(void)
sum = 0;
for (j = 0; j < count_unknowns; j++)
{
sum += inv_delta1[j] * array[i * max_column_count + j];
sum += inv_delta1[j] * my_array[i * max_column_count + j];
}
if (sum > array[(i * max_column_count) + count_unknowns] + toler)
if (sum > my_array[(i * max_column_count) + count_unknowns] + toler)
{
if (debug_inverse)
{
output_msg(sformatf(
"\tERROR: inequality not satisfied for %s, %e\n",
row_name[i],
(double) (sum - array[(i * max_column_count) + count_unknowns])));
(double) (sum - my_array[(i * max_column_count) + count_unknowns])));
}
rv = false;
}
@ -3475,7 +3475,7 @@ isotope_balance_equation(struct inverse *inv_ptr, int row, int n)
if (primary_jit == primary_ptr &&
jit->second.Get_isotope_number() == isotope_number)
{
array[row * max_column_count + i] +=
my_array[row * max_column_count + i] +=
f * jit->second.Get_total() * jit->second.Get_ratio();
}
}
@ -3500,7 +3500,7 @@ isotope_balance_equation(struct inverse *inv_ptr, int row, int n)
break;
}
column = col_epsilon + (k * inv_ptr->count_solns) + i;
array[row * max_column_count + column] +=
my_array[row * max_column_count + column] +=
f * jit->second.Get_ratio();
}
}
@ -3528,7 +3528,7 @@ isotope_balance_equation(struct inverse *inv_ptr, int row, int n)
(i * inv_ptr->count_isotope_unknowns) + k;
}
}
array[row * max_column_count + column] +=
my_array[row * max_column_count + column] +=
f * jit->second.Get_total();
}
}
@ -3548,11 +3548,11 @@ isotope_balance_equation(struct inverse *inv_ptr, int row, int n)
{
/* term for alpha phase unknowns */
column = col_phases + i;
array[row * max_column_count + column] =
my_array[row * max_column_count + column] =
isotope_ptr[j].ratio * isotope_ptr[j].coef;
/* term for phase isotope uncertainty unknown */
column = col_phase_isotopes + i * inv_ptr->count_isotopes + n;
array[row * max_column_count + column] = isotope_ptr[j].coef;
my_array[row * max_column_count + column] = isotope_ptr[j].coef;
break;
}
}
@ -3940,7 +3940,7 @@ phase_isotope_inequalities(struct inverse *inv_ptr)
{
for (k = 0; k < count_rows; k++)
{
array[k * max_column_count + column] = 0.0;
my_array[k * max_column_count + column] = 0.0;
}
continue;
}
@ -3948,7 +3948,7 @@ phase_isotope_inequalities(struct inverse *inv_ptr)
/*
* optimization
*/
array[(column - col_epsilon) * max_column_count + column] =
my_array[(column - col_epsilon) * max_column_count + column] =
SCALE_EPSILON /
inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
/*
@ -3957,17 +3957,17 @@ phase_isotope_inequalities(struct inverse *inv_ptr)
/* for phases constrained to precipitate */
if (inv_ptr->phases[i].constraint == PRECIPITATE)
{
array[count_rows * max_column_count + col_phases + i] =
my_array[count_rows * max_column_count + col_phases + i] =
inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
array[count_rows * max_column_count + column] = 1.0;
my_array[count_rows * max_column_count + column] = 1.0;
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
"iso pos");
row_name[count_rows] = string_hsave(token);
count_rows++;
array[count_rows * max_column_count + col_phases + i] =
my_array[count_rows * max_column_count + col_phases + i] =
inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + column] = -1.0;
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
"iso neg");
row_name[count_rows] = string_hsave(token);
@ -3977,17 +3977,17 @@ phase_isotope_inequalities(struct inverse *inv_ptr)
}
else if (inv_ptr->phases[i].constraint == DISSOLVE)
{
array[count_rows * max_column_count + col_phases + i] =
my_array[count_rows * max_column_count + col_phases + i] =
-inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
array[count_rows * max_column_count + column] = -1.0;
my_array[count_rows * max_column_count + column] = -1.0;
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
"iso pos");
row_name[count_rows] = string_hsave(token);
count_rows++;
array[count_rows * max_column_count + col_phases + i] =
my_array[count_rows * max_column_count + col_phases + i] =
-inv_ptr->phases[i].isotopes[j].ratio_uncertainty;
array[count_rows * max_column_count + column] = 1.0;
my_array[count_rows * max_column_count + column] = 1.0;
sprintf(token, "%s %s", inv_ptr->phases[i].phase->name,
"iso neg");
row_name[count_rows] = string_hsave(token);
@ -4460,13 +4460,13 @@ dump_netpath_pat(struct inverse *inv_ptr)
if (inv_ptr->pat == NULL)
return (OK);
array_save = array;
array_save = my_array;
l_delta_save = delta;
count_unknowns_save = count_unknowns;
max_row_count_save = max_row_count;
max_column_count_save = max_column_count;
array = NULL;
my_array = NULL;
delta = NULL;
count_unknowns = 0;
max_row_count = 0;
@ -4906,9 +4906,9 @@ dump_netpath_pat(struct inverse *inv_ptr)
count_pat_solutions);
}
//free_model_allocs();
array = (LDBLE *) free_check_null(array);
my_array = (LDBLE *) free_check_null(my_array);
delta = (LDBLE *) free_check_null(delta);
array = array_save;
my_array = array_save;
delta = l_delta_save;
count_unknowns = count_unknowns_save;
max_row_count = max_row_count_save;

1
src/phreeqcpp/kinetics.cpp
View File

@ -2450,6 +2450,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
* Set nsaver
*/
run_reactions_iterations = 0;
overall_iterations = 0;
kin_time_x = kin_time;
rate_kin_time = kin_time;
nsaver = i;

1
src/phreeqcpp/mainsubs.cpp
View File

@ -1085,6 +1085,7 @@ reactions(void)
rate_sim_time = 0;
for (reaction_step = 1; reaction_step <= count_steps; reaction_step++)
{
overall_iterations = 0;
sprintf(token, "Reaction step %d.", reaction_step);
if (reaction_step > 1 && incremental_reactions == FALSE)
{

207
src/phreeqcpp/model.cpp
View File

@ -96,6 +96,7 @@ model(void)
PhreeqcIWait(this);
#endif
iterations++;
overall_iterations++;
if (iterations > itmax - 1 && debug_model == FALSE
&& pr.logfile == TRUE)
{
@ -133,6 +134,7 @@ model(void)
{
jacobian_sums();
numerical_jacobian();
}
/*
* Full matrix with pure phases
@ -971,11 +973,11 @@ ineq(int in_kode)
{
for (j = 0; j < count_unknowns; j++)
{
array[j * (count_unknowns + 1) + i] = 0.0;
my_array[j * (count_unknowns + 1) + i] = 0.0;
}
for (j = 0; j < count_unknowns + 1; j++)
{
array[i * (count_unknowns + 1) + j] = 0.0;
my_array[i * (count_unknowns + 1) + j] = 0.0;
}
}
}
@ -1013,24 +1015,24 @@ ineq(int in_kode)
if (x[i]->type == SURFACE_CB1 && x[j]->type == SURFACE_CB2)
continue;
if (fabs(array[j * (count_unknowns + 1) + i]) > max)
if (fabs(my_array[j * (count_unknowns + 1) + i]) > max)
{
max = fabs(array[j * (count_unknowns + 1) + i]);
max = fabs(my_array[j * (count_unknowns + 1) + i]);
if (max > min_value)
break;
}
}
if (diagonal_scale == TRUE)
{
if (fabs(array[i * (count_unknowns + 1) + i]) < min_value)
if (fabs(my_array[i * (count_unknowns + 1) + i]) < min_value)
{
max = fabs(array[i * (count_unknowns + 1) + i]);
max = fabs(my_array[i * (count_unknowns + 1) + i]);
}
}
if (max == 0)
{
array[i * (count_unknowns + 1) + i] = 1e-5 * x[i]->moles;
my_array[i * (count_unknowns + 1) + i] = 1e-5 * x[i]->moles;
max = fabs(1e-5 * x[i]->moles);
}
}
@ -1041,11 +1043,11 @@ ineq(int in_kode)
min = 1e-12;
min = MIN_TOTAL;
array[x[i]->number * (count_unknowns + 1) + x[i]->number] += min;
my_array[x[i]->number * (count_unknowns + 1) + x[i]->number] += min;
if (fabs
(array[x[i]->number * (count_unknowns + 1) + x[i]->number]) <
(my_array[x[i]->number * (count_unknowns + 1) + x[i]->number]) <
min)
array[x[i]->number * (count_unknowns + 1) + x[i]->number] =
my_array[x[i]->number * (count_unknowns + 1) + x[i]->number] =
min;
max = 0.0;
@ -1059,9 +1061,9 @@ ineq(int in_kode)
x[j]->type != EXCH && x[j]->type != MH
&& x[j]->type != MH2O)
continue;
if (fabs(array[j * (count_unknowns + 1) + i]) > max)
if (fabs(my_array[j * (count_unknowns + 1) + i]) > max)
{
max = fabs(array[j * (count_unknowns + 1) + i]);
max = fabs(my_array[j * (count_unknowns + 1) + i]);
if (max > min_value)
break;
}
@ -1078,7 +1080,7 @@ ineq(int in_kode)
}
for (j = 0; j < count_unknowns; j++)
{
array[j * (count_unknowns + 1) + i] *= min_value / max;
my_array[j * (count_unknowns + 1) + i] *= min_value / max;
}
normal[i] = min_value / max;
}
@ -1151,7 +1153,7 @@ ineq(int in_kode)
{
/* Copy in saturation index equation (has mass or supersaturated) */
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
//if (it->second.Get_add_formula().size() == 0
@ -1185,7 +1187,7 @@ ineq(int in_kode)
* Alkalinity and solution phase boundary
*/
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
l_count_rows++;
@ -1196,7 +1198,7 @@ ineq(int in_kode)
else if (x[i]->type == GAS_MOLES && gas_in == TRUE)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
@ -1213,7 +1215,7 @@ ineq(int in_kode)
else if (x[i]->type == SS_MOLES && x[i]->ss_in == TRUE)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
res[l_count_rows] = 1.0;
@ -1258,6 +1260,7 @@ ineq(int in_kode)
continue;
if (x[i]->type == MH && pitzer_model == TRUE && pitzer_pe == FALSE)
continue;
if (mass_water_switch == TRUE && x[i] == mass_oxygen_unknown)
continue;
/*
@ -1294,7 +1297,7 @@ ineq(int in_kode)
}
}
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
if (mass_water_switch == TRUE && x[i] == mass_hydrogen_unknown)
@ -1303,7 +1306,7 @@ ineq(int in_kode)
for (j = 0; j < count_unknowns; j++)
{
ineq_array[l_count_rows * max_column_count + j] -=
2 * array[k * (count_unknowns + 1) + j];
2 * my_array[k * (count_unknowns + 1) + j];
}
}
l_count_rows++;
@ -1311,7 +1314,7 @@ ineq(int in_kode)
else if (x[i]->type == PITZER_GAMMA && full_pitzer == TRUE)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
l_count_rows++;
@ -1390,7 +1393,7 @@ ineq(int in_kode)
if (x[i]->type == MH2O)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
for (j = 0; j < count_unknowns; j++)
@ -1617,7 +1620,7 @@ ineq(int in_kode)
if (x[i]->type == MB && x[i]->moles < 0.0)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
(void *) &(my_array[i * (count_unknowns + 1)]),
(size_t) (count_unknowns + 1) * sizeof(LDBLE));
back_eq[l_count_rows] = i;
for (j = 0; j < count_unknowns; j++)
@ -1832,7 +1835,7 @@ ineq(int in_kode)
{
for (j = 0; j < count_unknowns; j++)
{
array[j * (count_unknowns + 1) + i] /= normal[i];
my_array[j * (count_unknowns + 1) + i] /= normal[i];
}
}
}
@ -1909,12 +1912,12 @@ jacobian_sums(void)
*/
for (i = 0; i < count_unknowns; i++)
{
array[i] = 0.0;
my_array[i] = 0.0;
}
for (i = 1; i < count_unknowns; i++)
{
memcpy((void *) &(array[i * (count_unknowns + 1)]),
(void *) &(array[0]), (size_t) count_unknowns * sizeof(LDBLE));
memcpy((void *) &(my_array[i * (count_unknowns + 1)]),
(void *) &(my_array[0]), (size_t) count_unknowns * sizeof(LDBLE));
}
/*
* Add constant terms
@ -1948,9 +1951,9 @@ jacobian_sums(void)
for (i = 0; i < count_unknowns; i++)
{
// using straight mu equation
array[mu_unknown->number * (count_unknowns + 1) + i] *= 0.5;
my_array[mu_unknown->number * (count_unknowns + 1) + i] *= 0.5;
}
array[mu_unknown->number * (count_unknowns + 1) +
my_array[mu_unknown->number * (count_unknowns + 1) +
mu_unknown->number] -= mass_water_aq_x;
}
/*
@ -1958,7 +1961,7 @@ jacobian_sums(void)
*/
if (mass_oxygen_unknown != NULL && mu_unknown != NULL)
{
array[mu_unknown->number * (count_unknowns + 1) +
my_array[mu_unknown->number * (count_unknowns + 1) +
mass_oxygen_unknown->number] -= mu_x * mass_water_aq_x;
}
/*
@ -1978,16 +1981,16 @@ jacobian_sums(void)
for (i = 0; i < count_unknowns; i++)
{
array[ah2o_unknown->number * (count_unknowns + 1) + i] *= factor;
my_array[ah2o_unknown->number * (count_unknowns + 1) + i] *= factor;
}
// activity of water term
array[ah2o_unknown->number * (count_unknowns + 1) +
my_array[ah2o_unknown->number * (count_unknowns + 1) +
ah2o_unknown->number] -= exp(s_h2o->la * LOG_10);
// mass of water term
if (mass_oxygen_unknown != NULL)
{
array[ah2o_unknown->number * (count_unknowns + 1) + mass_oxygen_unknown->number] -=
my_array[ah2o_unknown->number * (count_unknowns + 1) + mass_oxygen_unknown->number] -=
a*y_sum*(x_h2o*(0.5*tanh(lim) + 0.5) + (47.5*x_h2o - 50.0*a*y_sum)/(cosh(lim)*cosh(lim))) /
(x_h2o*x_h2o*x_h2o);
}
@ -1996,13 +1999,13 @@ jacobian_sums(void)
{
for (i = 0; i < count_unknowns; i++)
{
array[ah2o_unknown->number * (count_unknowns + 1) + i] *= -AH2O_FACTOR;
my_array[ah2o_unknown->number * (count_unknowns + 1) + i] *= -AH2O_FACTOR;
}
array[ah2o_unknown->number * (count_unknowns + 1) + ah2o_unknown->number] -=
my_array[ah2o_unknown->number * (count_unknowns + 1) + ah2o_unknown->number] -=
mass_water_aq_x * exp(s_h2o->la * LOG_10);
if (mass_oxygen_unknown != NULL)
{
array[ah2o_unknown->number * (count_unknowns + 1) + mass_oxygen_unknown->number] -=
my_array[ah2o_unknown->number * (count_unknowns + 1) + mass_oxygen_unknown->number] -=
(exp(s_h2o->la * LOG_10) - 1) * mass_water_aq_x;
}
}
@ -2019,7 +2022,7 @@ jacobian_sums(void)
//sqrt(8 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000) * tk_x *
// 1000);
sqrt(8 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000) * tk_x *
1000);
1000);
for (i = 0; i < count_unknowns; i++)
{
cxxSurfaceCharge *charge_ptr = NULL;
@ -2031,16 +2034,16 @@ jacobian_sums(void)
{
for (j = 0; j < count_unknowns; j++)
{
array[x[i]->number * (count_unknowns + 1) + j] *=
my_array[x[i]->number * (count_unknowns + 1) + j] *=
F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
charge_ptr->Get_grams());
}
array[x[i]->number * (count_unknowns + 1) + x[i]->number] -=
my_array[x[i]->number * (count_unknowns + 1) + x[i]->number] -=
sinh_constant * sqrt(mu_x) *
cosh(x[i]->master[0]->s->la * LOG_10);
if (mu_unknown != NULL)
{
array[x[i]->number * (count_unknowns + 1) +
my_array[x[i]->number * (count_unknowns + 1) +
mu_unknown->number] -=
0.5 * sinh_constant / sqrt(mu_x) *
sinh(x[i]->master[0]->s->la * LOG_10);
@ -2061,13 +2064,13 @@ jacobian_sums(void)
{
for (j = 0; j < count_unknowns; j++)
{
array[x[i]->number * (count_unknowns + 1) + j] *=
my_array[x[i]->number * (count_unknowns + 1) + j] *=
F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
charge_ptr->Get_grams());
}
array[x[i]->number * (count_unknowns + 1) + x[i]->number] -=
my_array[x[i]->number * (count_unknowns + 1) + x[i]->number] -=
charge_ptr->Get_capacitance0() * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ;
tk_x * LOG_10 / F_KJ_V_EQ;
}
}
}
@ -2322,6 +2325,12 @@ molalities(int allow_overflow)
}
if (dl_type_x != cxxSurface::NO_DL)
{
//s_h2o->tot_g_moles = s_h2o->moles;
//s_h2o->tot_dh2o_moles = 0.0;
if (calculating_deriv && use.Get_surface_ptr() != NULL) // DL_pitz
if (use.Get_surface_ptr()->Get_debye_lengths() > 0)
s_h2o->moles = mass_water_bulk_x / gfw_water;
s_h2o->tot_g_moles = s_h2o->moles;
s_h2o->tot_dh2o_moles = 0.0;
}
@ -2377,8 +2386,9 @@ molalities(int allow_overflow)
* other terms for diffuse layer model
*/
if (use.Get_surface_ptr() != NULL
&& use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC
&& dl_type_x != cxxSurface::NO_DL)
&& dl_type_x != cxxSurface::NO_DL
&& (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC
|| pitzer_model || calculating_deriv)) // DL_pitz
{
calc_all_donnan();
}
@ -3724,7 +3734,9 @@ reset(void)
else if (x[i]->type == MH2O)
{
if (mass_water_switch == TRUE)
{
continue;
}
/*if (fabs(delta[i]) > epsilon * mass_water_aq_x) converge=FALSE; */
/* ln(gh2o) + delta, log(gh2o) + d, gh2o * 10**d */
d = exp(delta[i]);
@ -3738,7 +3750,6 @@ reset(void)
(double) delta[i], "10**d/c", (double) d));
}
mass_water_aq_x *= d;
mass_water_bulk_x = mass_water_aq_x + mass_water_surfaces_x;
if (debug_model == TRUE && dl_type_x != cxxSurface::NO_DL)
{
@ -4293,27 +4304,22 @@ residuals(void)
}
else
{
/*
* sinh_constant is (8 e e0 R T 1000)**1/2
* = sqrt(8*EPSILON*EPSILON_ZERO*(R_KJ_DEG_MOL*1000)*t_x*1000)
* ~ 0.1174 at 25C
*/
residual[i] =
sinh_constant * sqrt(mu_x) *
sinh(x[i]->master[0]->s->la * LOG_10) -
x[i]->f * F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
residual[i] = sinh_constant * sqrt(mu_x) * sinh(x[i]->master[0]->s->la * LOG_10) -
x[i]->f * F_C_MOL / (charge_ptr->Get_specific_area() * charge_ptr->Get_grams());
}
if (debug_model == TRUE)
{
output_msg(sformatf( "Charge/Potential\n"));
if (charge_ptr->Get_grams() > 0)
{
output_msg(sformatf(
"\tSum of surface charge %e eq\n",
(double) (x[i]->f
/* F_C_MOL / (x[i]->surface_charge->specific_area * x[i]->surface_charge->grams) */
)));
if (dl_type_x != cxxSurface::NO_DL)
output_msg(sformatf("\tSum of surface + diffuse layer charge %e eq\n", (double)(x[i]->f)));
else
output_msg(sformatf(
"\tSum of surface charge %e eq\n",
(double)(x[i]->f
/* F_C_MOL / (x[i]->surface_charge->specific_area * x[i]->surface_charge->grams) */
)));
}
else
{
@ -4432,7 +4438,10 @@ residuals(void)
output_msg(sformatf( "Charge/Potential\n"));
if (charge_ptr->Get_grams() > 0)
{
output_msg(sformatf(
if (dl_type_x != cxxSurface::NO_DL)
output_msg(sformatf("\tSum of surface + diffuse layer charge %e eq\n", (double)(x[i]->f)));
else
output_msg(sformatf(
"\tSum of surface charge %e eq\n",
(double) (x[i]->f
/* F_C_MOL / (x[i]->surface_charge->specific_area * x[i]->surface_charge->grams) */
@ -4669,7 +4678,7 @@ residuals(void)
/*
* Store residuals in array
*/
array[(i + 1) * (count_unknowns + 1) - 1] = residual[i];
my_array[(i + 1) * (count_unknowns + 1) - 1] = residual[i];
sum_residual += fabs(residual[i]);
}
/*
@ -4800,9 +4809,9 @@ initial_guesses(void)
}
}
else if (x[i]->type == SURFACE_CB)
{
x[i]->master[0]->s->la = 0.0;
}
{
x[i]->master[0]->s->la = 0.0;
}
}
return (OK);
}
@ -5216,6 +5225,7 @@ surface_model(void)
}
if (model() == ERROR)
return (ERROR);
g_iterations = 0;
if (use.Get_surface_ptr()->Get_dl_type() == cxxSurface::DONNAN_DL)
@ -5230,6 +5240,7 @@ surface_model(void)
mb_sums();
if (model() == ERROR)
return (ERROR);
if (!use.Get_surface_ptr()->Get_related_phases()
&& !use.Get_surface_ptr()->Get_related_rate())
initial_surface_water();
@ -5303,7 +5314,7 @@ free_model_allocs(void)
}
x = (struct unknown **) free_check_null(x);
max_unknowns = 0;
array = (LDBLE *) free_check_null(array);
my_array = (LDBLE *) free_check_null(my_array);
delta = (LDBLE *) free_check_null(delta);
residual = (LDBLE *) free_check_null(residual);
s_x = (struct species **) free_check_null(s_x);
@ -5484,12 +5495,12 @@ numerical_jacobian(void)
for (i = 0; i < count_unknowns; i++)
{
array[i] = 0.0;
my_array[i] = 0.0;
}
for (i = 1; i < count_unknowns; i++)
{
memcpy((void *) &(array[i * (count_unknowns + 1)]),
(void *) &(array[0]), (size_t) count_unknowns * sizeof(LDBLE));
memcpy((void *) &(my_array[i * (count_unknowns + 1)]),
(void *) &(my_array[0]), (size_t) count_unknowns * sizeof(LDBLE));
}
base = (LDBLE *) PHRQ_malloc((size_t) count_unknowns * sizeof(LDBLE));
@ -5500,7 +5511,7 @@ numerical_jacobian(void)
base[i] = residual[i];
}
d = 0.0001;
d1 = d * log(10.0);
d1 = d * LOG_10;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
{
@ -5516,24 +5527,32 @@ numerical_jacobian(void)
case SURFACE_CB1:
case SURFACE_CB2:
x[i]->master[0]->s->la += d;
d2 = d1;
d2 = d * LOG_10;
break;
case MH:
s_eminus->la += d;
d2 = d1;
d2 = d * LOG_10;
break;
case AH2O:
x[i]->master[0]->s->la += d;
d2 = d1;
d2 = d * LOG_10;
break;
case PITZER_GAMMA:
x[i]->s->lg += d;
d2 = d;
break;
case MH2O:
mass_water_aq_x *= (1.0 + d);
//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 = log(1.0 + d);
d2 = d1;
break;
case MU:
d2 = d * mu_x;
@ -5558,8 +5577,8 @@ numerical_jacobian(void)
delta[j] = 0.0;
}
/*d2 = -1e-8; */
d2 = -d * x[i]->moles;
d2 = -.1 * x[i]->moles;
d2 = d * 10 * x[i]->moles;
//d2 = -.1 * x[i]->moles;
/*
if (d2 > -1e-10) d2 = -1e-10;
calculating_deriv = FALSE;
@ -5596,11 +5615,13 @@ numerical_jacobian(void)
LDBLE t = (LDBLE) pow((LDBLE) 10.0, (LDBLE) (DBL_MAX_10_EXP - 50.0));
if (residual[j] > t)
{
array[j * (count_unknowns + 1) + i] = -pow(10.0, DBL_MAX_10_EXP - 50.0);
my_array[j * (count_unknowns + 1) + i] = -pow(10.0, DBL_MAX_10_EXP - 50.0);
}
else
{
array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
my_array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
if (x[i]->type == MH2O) // DL_pitz
my_array[j * (count_unknowns + 1) + i] *= mass_water_aq_x;
}
}
else if (residual[j] < -1.0e101)
@ -5608,17 +5629,21 @@ numerical_jacobian(void)
LDBLE t = pow((LDBLE) 10.0, (LDBLE) (DBL_MIN_10_EXP + 50.0));
if (residual[j] < -t)
{
array[j * (count_unknowns + 1) + i] = pow(10.0, DBL_MIN_10_EXP + 50.0);
my_array[j * (count_unknowns + 1) + i] = pow(10.0, DBL_MIN_10_EXP + 50.0);
}
else
{
array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
my_array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
if (x[i]->type == MH2O) // DL_pitz
my_array[j * (count_unknowns + 1) + i] *= mass_water_aq_x;
}
}
else
{
array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
if (!PHR_ISFINITE(array[j * (count_unknowns + 1) + i]))
my_array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
if (x[i]->type == MH2O) // DL_pitz
my_array[j * (count_unknowns + 1) + i] *= mass_water_aq_x;
if (!PHR_ISFINITE(my_array[j * (count_unknowns + 1) + i]))
{
//fprintf(stderr, "oops, got NaN: %e, %e, %e, %e\n", residual[j], base[j], d2, array[j * (count_unknowns + 1) + i]);
}
@ -5642,10 +5667,10 @@ numerical_jacobian(void)
break;
case MH:
s_eminus->la -= d;
if (array[i * (count_unknowns + 1) + i] == 0)
if (my_array[i * (count_unknowns + 1) + i] == 0)
{
/*output_msg(sformatf( "Zero diagonal for MH\n")); */
array[i * (count_unknowns + 1) + i] =
my_array[i * (count_unknowns + 1) + i] =
under(s_h2->lm) * 2;
}
break;
@ -5653,7 +5678,13 @@ numerical_jacobian(void)
x[i]->s->lg -= d;
break;
case MH2O:
mass_water_aq_x /= (1 + d);
//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 -= 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;
break;
case MU:

105
src/phreeqcpp/pitzer.cpp
View File

@ -908,7 +908,7 @@ pitzer(void)
C INITIALIZE
C
*/
CONV = 1.0 / log(10.0);
CONV = 1.0 / LOG_10;
XX = 0.0;
OSUM = 0.0;
/*n
@ -1243,7 +1243,7 @@ pitzer(void)
C INITIALIZE
C
*/
CONV = 1.0 / log(10.0);
CONV = 1.0 / LOG_10;
XX = 0.0;
OSUM = 0.0;
I = mu_x;
@ -1731,6 +1731,7 @@ set_pz(int initial)
tk_x = tc_x + 273.15;
patm_x = solution_ptr->Get_patm(); // done in calc_rho_0(tc, pa)
potV_x = solution_ptr->Get_potV(); // added in DL_pitz
/*
* H+, e-, H2O
@ -1834,7 +1835,8 @@ pitzer_revise_guesses(void)
LDBLE weight, f;
max_iter = 100;
/* gammas(mu_x); */
if (iterations < 0 && (use.Get_surface_in() || use.Get_exchange_in()))
gammas_pz(); // DL_pitz : for SURF estimates
l_iter = 0;
repeat = TRUE;
fail = FALSE;
@ -1916,7 +1918,7 @@ pitzer_revise_guesses(void)
{
repeat = TRUE;
x[i]->master[0]->s->la += logd;
/*!!!!*/ if (x[i]->master[0]->s->la < -999.)
if (x[i]->master[0]->s->la < -999.)
x[i]->master[0]->s->la = MIN_RELATED_LOG_ACTIVITY;
}
else if (f > d * fabs(x[i]->moles)
@ -1975,7 +1977,7 @@ pitzer_revise_guesses(void)
{
mu_x = 1e-8;
}
/*gammas(mu_x); */
//gammas_pz();
return (OK);
}
@ -1983,11 +1985,12 @@ pitzer_revise_guesses(void)
int Phreeqc::
jacobian_pz(void)
/* ---------------------------------------------------------------------- */
{
{ // calculate the derivatives numerically
LDBLE *base;
LDBLE d, d1, d2;
int i, j;
calculating_deriv = 1;
Restart:
int pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x);
@ -2005,7 +2008,7 @@ Restart:
base[i] = residual[i];
}
d = 0.0001;
d1 = d * log(10.0);
d1 = d * LOG_10;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
{
@ -2021,11 +2024,13 @@ Restart:
case SURFACE_CB1:
case SURFACE_CB2:
x[i]->master[0]->s->la += d;
d2 = d1;
//d2 = d1;
d2 = d * LOG_10;
break;
case AH2O:
x[i]->master[0]->s->la += d;
d2 = d1;
//d2 = d1;
d2 = d * LOG_10;
break;
case PITZER_GAMMA:
if (!full_pitzer)
@ -2034,15 +2039,25 @@ Restart:
d2 = d;
break;
case MH2O:
mass_water_aq_x *= (1.0 + d);
//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 = log(1.0 + d);
//d2 = log(1.0 + d);
d2 = d1;
break;
case MH:
if (pitzer_pe == TRUE)
{
s_eminus->la += d;
d2 = d1;
//d2 = d1;
d2 = d * LOG_10;
break;
}
else
@ -2065,9 +2080,22 @@ Restart:
gammas(mu_x);
break;
case PP:
case SS_MOLES:
continue;
break;
case SS_MOLES:
//continue;
//break;
if (x[i]->ss_in == FALSE)
continue;
for (j = 0; j < count_unknowns; j++)
{
delta[j] = 0.0;
}
d2 = d * 10 * x[i]->moles;
delta[i] = d2;
reset();
d2 = delta[i];
break;
}
molalities(TRUE);
if (max_unknowns > pz_max_unknowns)
@ -2083,8 +2111,9 @@ Restart:
residuals();
for (j = 0; j < count_unknowns; j++)
{
array[j * (count_unknowns + 1) + i] =
-(residual[j] - base[j]) / d2;
my_array[j * (count_unknowns + 1) + i] = -(residual[j] - base[j]) / d2;
if (x[i]->type == MH2O) // DL_pitz
my_array[j * (count_unknowns + 1) + i] *= mass_water_aq_x;
}
switch (x[i]->type)
{
@ -2102,9 +2131,9 @@ Restart:
break;
case MH:
s_eminus->la -= d;
if (array[i * (count_unknowns + 1) + i] == 0)
if (my_array[i * (count_unknowns + 1) + i] == 0)
{
array[i * (count_unknowns + 1) + i] =
my_array[i * (count_unknowns + 1) + i] =
exp(s_h2->lm * LOG_10) * 2;
}
break;
@ -2112,7 +2141,13 @@ Restart:
x[i]->s->lg -= d;
break;
case MH2O:
mass_water_aq_x /= (1 + d);
//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 -= 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;
break;
case MU:
@ -2125,7 +2160,10 @@ Restart:
continue;
x[i]->moles -= d2;
break;
case SS_MOLES:
delta[i] = -d2;
reset();
break;
}
}
molalities(TRUE);
@ -2134,6 +2172,7 @@ Restart:
mb_sums();
residuals();
free_check_null(base);
calculating_deriv = 0;
return OK;
}
@ -2217,6 +2256,7 @@ model_pz(void)
PhreeqcIWait(this);
#endif
iterations++;
overall_iterations++;
if (iterations > itmax - 1 && debug_model == FALSE
&& pr.logfile == TRUE)
{
@ -2307,13 +2347,13 @@ model_pz(void)
reprep();
full_pitzer = false;
}
/* debug
species_list_sort();
sum_species();
print_species();
print_exchange();
print_surface();
*/
//debug
//species_list_sort();
//sum_species();
//print_species();
//print_exchange();
//print_surface();
if (stop_program == TRUE)
{
break;
@ -2422,7 +2462,7 @@ gammas_pz()
* Need exchange gammas for pitzer
*/
int i, j;
LDBLE coef;
LDBLE coef, equiv;
/* Initialize */
k_temp(tc_x, patm_x);
/*
@ -2456,9 +2496,18 @@ gammas_pz()
break;
}
}
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) // DL_pitz
{
/* mole fraction */
equiv = 1.0;
}
else
{
equiv = s_x[i]->equiv;
}
if (s_x[i]->alk > 0)
{
s_x[i]->lg = log10(s_x[i]->equiv / s_x[i]->alk);
s_x[i]->lg = log10(equiv / s_x[i]->alk);
s_x[i]->dg = 0.0;
}
else

136
src/phreeqcpp/prep.cpp
View File

@ -89,10 +89,10 @@ prep(void)
residual = (LDBLE *) PHRQ_malloc( (size_t) count_unknowns * sizeof( LDBLE ));
if (residual == NULL) malloc_error();
*/
array =
my_array =
(LDBLE *) PHRQ_malloc((size_t) (max_unknowns + 1) *
max_unknowns * sizeof(LDBLE));
if (array == NULL)
if (my_array == NULL)
malloc_error();
delta = (LDBLE *) PHRQ_malloc((size_t) max_unknowns * sizeof(LDBLE));
if (delta == NULL)
@ -120,12 +120,26 @@ prep(void)
*/
quick_setup();
}
if (debug_mass_balance)
{
output_msg(sformatf("\nTotals for the equation solver.\n"));
output_msg(sformatf("\n\tRow\tName Type Total moles\n"));
for (int i = 0; i < count_unknowns; i++)
{
if (x[i]->type == PITZER_GAMMA)
continue;
output_msg(sformatf("\t%3d\t%-17s%2d %15.6e\n",
x[i]->number, x[i]->description, (int)x[i]->type, (double)x[i]->moles));
}
output_msg(sformatf("\n\n"));
}
if (get_input_errors() > 0)
{
error_msg("Program stopping due to input errors.", STOP);
}
if (sit_model) sit_make_lists();
if (pitzer_model) pitzer_make_lists();
if (pitzer_model)
pitzer_make_lists();
return (OK);
}
@ -420,7 +434,7 @@ build_gas_phase(void)
*/
if (debug_prep == TRUE)
{
output_msg(sformatf( "\n\tMass balance summations %s.\n\n",
output_msg(sformatf( "\n\tMass balance summations. %s.\n",
phase_ptr->name));
}
@ -543,27 +557,27 @@ build_gas_phase(void)
}
col = master_ptr->unknown->number;
coef = coef_elt * rxn_ptr->coef;
store_jacob(&(phase_ptr->moles_x),
&(array[row + col]), coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
master_ptr->s->name, (double) coef,
row / (count_unknowns + 1), col));
}
store_jacob(&(phase_ptr->moles_x),
&(my_array[row + col]), coef);
}
if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE)
{
/* derivative wrt total moles of gas */
store_jacob(&(phase_ptr->fraction_x),
&(array[row + gas_unknown->number]), coef_elt);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
"gas moles", (double) elt_list[j].coef,
row / (count_unknowns + 1),
gas_unknown->number));
}
store_jacob(&(phase_ptr->fraction_x),
&(my_array[row + gas_unknown->number]), coef_elt);
}
}
/*
@ -634,13 +648,13 @@ build_gas_phase(void)
}
col = master_ptr->unknown->number;
coef = rxn_ptr->coef;
store_jacob(&(phase_ptr->p_soln_x), &(array[row + col]), coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
master_ptr->s->name, (double) coef,
row / (count_unknowns + 1), col));
}
store_jacob(&(phase_ptr->p_soln_x), &(my_array[row + col]), coef);
}
}
}
@ -733,11 +747,11 @@ build_ss_assemblage(void)
{
col = x[i]->number - 1;
}
store_jacob(&(x[i]->phase->dnc), &(array[row + col]), -1);
store_jacob(&(x[i]->phase->dnc), &(my_array[row + col]), -1);
/* next dnb terms */
col++;
store_jacob(&(x[i]->phase->dnb), &(array[row + col]), -1);
store_jacob(&(x[i]->phase->dnb), &(my_array[row + col]), -1);
}
else
{
@ -750,12 +764,12 @@ build_ss_assemblage(void)
if ((int) j != x[i]->ss_comp_number)
{
/* store_jacob (&(s_s_ptr->dn), &(array[row + col + j]), -1.0); */
store_jacob(&(x[i]->phase->dn), &(array[row + col + j]),
store_jacob(&(x[i]->phase->dn), &(my_array[row + col + j]),
-1.0);
}
else
{
store_jacob(&(x[i]->phase->dnb), &(array[row + col + j]),
store_jacob(&(x[i]->phase->dnb), &(my_array[row + col + j]),
-1.0);
}
}
@ -892,7 +906,7 @@ build_jacobian_sums(int k)
}
coef = mb_unknowns[i].coef;
if (debug_prep == TRUE)
output_msg(sformatf( "\n\tMass balance eq: %s\t%f\n",
output_msg(sformatf( "\n\tMass balance eq: %-13s\t%f\trow\tcol\n",
mb_unknowns[i].unknown->description, (double) coef));
store_dn(k, mb_unknowns[i].source, mb_unknowns[i].unknown->number,
coef, mb_unknowns[i].gamma_source);
@ -912,17 +926,17 @@ build_jacobian_sums(int k)
/* term for water, sum of all surfaces */
source = &s[k]->tot_dh2o_moles;
target =
&(array
&(my_array
[mb_unknowns[i].unknown->number * (count_unknowns + 1) +
mass_oxygen_unknown->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
"sum[dn(i,s)/dlnwater]", (double) coef,
mb_unknowns[i].unknown->number,
mass_oxygen_unknown->number));
}
store_jacob(source, target, coef);
}
/* terms for psi, one for each surface */
@ -933,15 +947,15 @@ build_jacobian_sums(int k)
continue;
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
source = s_diff_layer[k][charge_ptr->Get_name()].Get_dx_moles_address();
target = &(array[mb_unknowns[i].unknown->number *
target = &(my_array[mb_unknowns[i].unknown->number *
(count_unknowns + 1) + x[j]->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
"dg/dlny", (double) coef,
mb_unknowns[i].unknown->number, x[j]->number));
}
store_jacob(source, target, coef);
count_g++;
if (count_g >= (int) use.Get_surface_ptr()->Get_surface_charges().size())
break;
@ -972,17 +986,17 @@ build_jacobian_sums(int k)
if (kk >= 0)
{
source = s_diff_layer[k][charge_ptr->Get_name()].Get_drelated_moles_address();
target = &(array[mb_unknowns[i].unknown->number *
target = &(my_array[mb_unknowns[i].unknown->number *
(count_unknowns + 1) + x[kk]->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf(
"\t\t%-24s%10.3f\t%d\t%d\n", "dphase",
"\t\t%-24s%10.3f\t%d\t%d", "dphase",
(double) coef,
mb_unknowns[i].unknown->number,
x[kk]->number));
}
store_jacob(source, target, coef);
}
count_g++;
if (count_g >= (int) use.Get_surface_ptr()->Get_surface_charges().size())
@ -1001,17 +1015,17 @@ build_jacobian_sums(int k)
if (mb_unknowns[i].unknown->number == x[j]->number)
{
source = s_diff_layer[k][charge_ptr->Get_name()].Get_dx_moles_address();
target = &(array[mb_unknowns[i].unknown->number *
target = &(my_array[mb_unknowns[i].unknown->number *
(count_unknowns + 1) + x[j]->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf(
"\t\t%-24s%10.3f\t%d\t%d\n", "dg/dlny",
"\t\t%-24s%10.3f\t%d\t%d", "dg/dlny",
(double) coef,
mb_unknowns[i].unknown->number,
x[j]->number));
}
store_jacob(source, target, coef);
/* term for related phase */
/* has related phase */
@ -1030,17 +1044,17 @@ build_jacobian_sums(int k)
if (kk >= 0)
{
source = s_diff_layer[k][charge_ptr->Get_name()].Get_drelated_moles_address();
target = &(array[mb_unknowns[i].unknown->number *
target = &(my_array[mb_unknowns[i].unknown->number *
(count_unknowns + 1) + x[kk]->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf(
"\t\t%-24s%10.3f\t%d\t%d\n",
"\t\t%-24s%10.3f\t%d\t%d",
"dphase", (double) coef,
mb_unknowns[i].unknown->number,
x[kk]->number));
}
store_jacob(source, target, coef);
}
}
@ -1048,18 +1062,18 @@ build_jacobian_sums(int k)
{
/* term for water, for same surfaces */
source = s_diff_layer[k][charge_ptr->Get_name()].Get_dh2o_moles_address();
target = &(array[mb_unknowns[i].unknown->number *
target = &(my_array[mb_unknowns[i].unknown->number *
(count_unknowns + 1) +
mass_oxygen_unknown->number]);
store_jacob(source, target, coef);
if (debug_prep == TRUE)
{
output_msg(sformatf(
"\t\t%-24s%10.3f\t%d\t%d\n",
"\t\t%-24s%10.3f\t%d\t%d",
"dn(i,s)/dlnwater", (double) coef,
mb_unknowns[i].unknown->number,
mass_oxygen_unknown->number));
}
store_jacob(source, target, coef);
}
break;
}
@ -1102,7 +1116,7 @@ build_mb_sums(void)
if (debug_prep == TRUE)
{
output_msg(sformatf( "\n\tMass balance summations.\n\n"));
output_msg(sformatf( "\n\tMass balance summations.\n"));
}
for (i = 0; i < count_mb_unknowns; i++)
{
@ -1236,7 +1250,7 @@ build_model(void)
{
s[i]->dz[j] = s[i]->rxn_x->dz[j];
}
if (debug_prep == TRUE)
if (debug_mass_action == TRUE)
{
output_msg(sformatf( "\n%s\n\tMass-action equation\n",
s[i]->name));
@ -1265,7 +1279,7 @@ build_model(void)
}
if (debug_prep == TRUE)
{
output_msg(sformatf( "\tElement composition %s\n",
output_msg(sformatf( "\n%s, Element composition:\n",
trxn.token[0].s->name));
for (j = 0; j < count_elts; j++)
{
@ -1274,12 +1288,12 @@ build_model(void)
(double) elt_list[j].coef));
}
}
if (debug_prep == TRUE)
{
output_msg(sformatf( "\n\tMass balance equation\n",
s[i]->name));
trxn_print();
}
//if (debug_prep == TRUE)
//{
// output_msg(sformatf( "\n\tMass balance equation\n",
// s[i]->name));
// trxn_print();
//}
if (s[i]->type < EMINUS)
{
mb_for_species_aq(i);
@ -1300,6 +1314,7 @@ build_model(void)
build_mb_sums();
}
#endif
if (!pitzer_model && !sit_model)
build_jacobian_sums(i);
/*
@ -1317,7 +1332,7 @@ build_model(void)
build_species_list(i);
}
}
if (dl_type_x != cxxSurface::NO_DL && (pitzer_model == TRUE || sit_model == TRUE))
if (dl_type_x != cxxSurface::NO_DL && (/*pitzer_model == TRUE || */sit_model == TRUE)) //DL_pitz
{
error_msg("-diffuse_layer option not available for Pizer or SIT model",
STOP);
@ -1344,9 +1359,8 @@ build_model(void)
}
}
/*
* For Pizer model add lg unknown for each aqueous species
* For Pitzer model add lg unknown for each aqueous species
*/
if (pitzer_model == TRUE || sit_model == TRUE)
{
j0 = count_unknowns;
@ -1366,7 +1380,7 @@ build_model(void)
count_unknowns++;
}
}
/*
/*
* Rewrite phases to current master species
*/
for (i = 0; i < count_phases; i++)
@ -4880,14 +4894,14 @@ store_dn(int k, LDBLE * source, int row, LDBLE coef_in, LDBLE * gamma_source)
{
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
"Activity coefficient", (double) (-1.0 * coef_in),
row / (count_unknowns + 1), mu_unknown->number));
}
/* mu term */
if (gamma_source != NULL)
{
store_jacob(gamma_source, &array[row + mu_unknown->number],
store_jacob(gamma_source, &my_array[row + mu_unknown->number],
-1.0 * coef_in);
}
}
@ -4898,12 +4912,12 @@ store_dn(int k, LDBLE * source, int row, LDBLE coef_in, LDBLE * gamma_source)
{
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
mass_oxygen_unknown->master[0]->s->name,
(double) coef_in, row / (count_unknowns + 1),
mass_oxygen_unknown->number));
}
store_jacob(source, &(array[row + mass_oxygen_unknown->number]),
store_jacob(source, &(my_array[row + mass_oxygen_unknown->number]),
coef_in);
}
if (s[k] == s_h2o)
@ -4919,21 +4933,21 @@ store_dn(int k, LDBLE * source, int row, LDBLE coef_in, LDBLE * gamma_source)
{
master_ptr = rxn_ptr->s->primary;
}
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%s\n", master_ptr->s->name));
}
//if (debug_prep == TRUE)
//{
// output_msg(sformatf( "\t\t%s\n", master_ptr->s->name));
//}
if (master_ptr == NULL ||master_ptr->unknown == NULL)
continue;
col = master_ptr->unknown->number;
coef = coef_in * rxn_ptr->coef;
store_jacob(source, &(array[row + col]), coef);
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d\n",
output_msg(sformatf( "\t\t%-24s%10.3f\t%d\t%d",
master_ptr->s->name, (double) coef,
row / (count_unknowns + 1), col));
}
store_jacob(source, &(my_array[row + col]), coef);
}
return (OK);
}
@ -4952,7 +4966,7 @@ store_jacob(LDBLE * source, LDBLE * target, LDBLE coef)
{
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\tjacob1 %d\n", count_sum_jacob1));
output_msg(sformatf( "\tjacob1 %d\n", count_sum_jacob1));
}
sum_jacob1[count_sum_jacob1].source = source;
sum_jacob1[count_sum_jacob1++].target = target;
@ -4967,7 +4981,7 @@ store_jacob(LDBLE * source, LDBLE * target, LDBLE coef)
{
if (debug_prep == TRUE)
{
output_msg(sformatf( "\t\tjacob2 %d\n", count_sum_jacob2));
output_msg(sformatf( "\tjacob2 %d\n", count_sum_jacob2));
}
sum_jacob2[count_sum_jacob2].source = source;
sum_jacob2[count_sum_jacob2].target = target;
@ -4991,7 +5005,7 @@ store_jacob0(int row, int column, LDBLE coef)
* Stores in list a constant coef which will be added into jacobian array
*/
sum_jacob0[count_sum_jacob0].target =
&(array[row * (count_unknowns + 1) + column]);
&(my_array[row * (count_unknowns + 1) + column]);
sum_jacob0[count_sum_jacob0++].coef = coef;
/* Check space */
if (count_sum_jacob0 >= max_sum_jacob0)

10
src/phreeqcpp/print.cpp
View File

@ -195,7 +195,8 @@ punch_all(void)
/*
* new line for punch_file
*/
punch_msg("\n");
if (current_selected_output->Get_new_line())
punch_msg("\n");
/*
* signal end of row
@ -284,7 +285,7 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
output_msg(sformatf(
"\n\tSpecies \t Moles \tMoles excess\t g\n"));
}
if ((mass_water_surface = charge_ptr->Get_mass_water()))
if (mass_water_surface = charge_ptr->Get_mass_water())
{
count_elts = 0;
paren_count = 0;
@ -2320,7 +2321,10 @@ print_totals(void)
output_msg(sformatf("%45s%6.2f\n",
"Percent error, 100*(Cat-|An|)/(Cat+|An|) = ",
(double) (100 * cb_x / total_ions_x)));
output_msg(sformatf("%45s%3d\n", "Iterations = ", iterations));
if (iterations == overall_iterations)
output_msg(sformatf("%45s%3d\n", "Iterations = ", iterations));
else
output_msg(sformatf("%45s%3d (%d overall)\n", "Iterations = ", iterations, overall_iterations));
if (pitzer_model == TRUE || sit_model == TRUE)
{
if (always_full_pitzer == FALSE)

24
src/phreeqcpp/read.cpp
View File

@ -4749,9 +4749,10 @@ read_selected_output(void)
"isotopes", /* 46 */
"calculate_values", /* 47 */
"equilibrium_phase", /* 48 */
"active" /* 49 */
"active", /* 49 */
"new_line" /* 50 */
};
int count_opt_list = 50;
int count_opt_list = 51;
int i, l;
char file_name[MAX_LENGTH], token[MAX_LENGTH];
@ -4777,6 +4778,7 @@ read_selected_output(void)
// n_user = 1, old definition, keep old definition
SelectedOutput & so_ref = so->second;
temp_selected_output.Set_active ( so_ref.Get_active() );
temp_selected_output.Set_new_line ( so_ref.Get_new_line() );
temp_selected_output.Set_inverse ( so_ref.Get_inverse() );
temp_selected_output.Set_sim ( so_ref.Get_sim() );
temp_selected_output.Set_state ( so_ref.Get_state() );
@ -5129,6 +5131,12 @@ read_selected_output(void)
}
opt_save = OPTION_ERROR;
break;
case 50: /* new_line */
temp_selected_output.Set_new_def(true);
value = get_true_false(next_char, TRUE);
temp_selected_output.Set_new_line(value != FALSE);
opt_save = OPTION_ERROR;
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
@ -8367,9 +8375,11 @@ read_debug(void)
"force_numerical_fixed_volume", /* 19 */
"equi_delay", /* 20 */
"minimum_total", /* 21 */
"min_total" /* 22 */
"min_total", /* 22 */
"debug_mass_action", /* 23 */
"debug_mass_balance" /* 24 */
};
int count_opt_list = 23;
int count_opt_list = 25;
/*
* Read parameters:
* ineq_tol;
@ -8471,6 +8481,12 @@ read_debug(void)
MIN_TOTAL_SS = MIN_TOTAL/100;
MIN_RELATED_SURFACE = MIN_TOTAL*100;
break;
case 23: /* debug_mass_action */
debug_mass_action = get_true_false(next_char, TRUE);
break;
case 24: /* debug_mass_balance */
debug_mass_balance = get_true_false(next_char, TRUE);
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;

10
src/phreeqcpp/readtr.cpp
View File

@ -817,7 +817,7 @@ read_transport(void)
*/
if (count_por == 0)
{
if (old_cells < all_cells && multi_Dflag && simul_tr == 1)
if (old_cells < all_cells && multi_Dflag /*&& simul_tr == 1*/)
{
multi_Dpor = (multi_Dpor < 1e-10 ? 1e-10 : multi_Dpor);
if (multi_Dpor > 1e-10)
@ -827,7 +827,11 @@ read_transport(void)
error_string = sformatf(
"No porosities were read; set to minimal value of 1e-10 for -multi_D.");
warning_msg(error_string);
for (i = 1; i < all_cells; i++)
if (simul_tr == 1)
j = 1;
else
j = old_cells + 1;
for (i = j; i < all_cells; i++)
cell_data[i].por = multi_Dpor;
}
}
@ -847,7 +851,7 @@ read_transport(void)
{
for (i = 1; i <= count_por; i++)
cell_data[i].por = pors[i - 1];
if (max_cells > count_por)
if (all_cells > count_por)
{
error_string = sformatf(
"Porosities were read for %d cells. Last value is used till cell %d.",

21
src/phreeqcpp/sit.cpp
View File

@ -312,7 +312,7 @@ sit(void)
C INITIALIZE
C
*/
//CONV = 1.0 / log(10.0);
//CONV = 1.0 / LOG_10;
XI = 0.0e0;
XX = 0.0e0;
OSUM = 0.0e0;
@ -367,7 +367,7 @@ sit(void)
C
*/
AGAMMA = 3*sit_A0; /* Grenthe p 379 */
A = AGAMMA / log(10.0);
A = AGAMMA / LOG_10;
/*
* F is now for log10 gamma
*/
@ -445,7 +445,7 @@ sit(void)
C
*/
/*COSMOT = 1.0e0 + 2.0e0 * OSMOT / OSUM;*/
COSMOT = 1.0e0 + OSMOT*log(10.0) / OSUM;
COSMOT = 1.0e0 + OSMOT*LOG_10 / OSUM;
/*
C
C CALCULATE THE ACTIVITY OF WATER
@ -490,7 +490,7 @@ sit(void)
C INITIALIZE
C
*/
//CONV = 1.0 / log(10.0);
//CONV = 1.0 / LOG_10;
XI = 0.0e0;
XX = 0.0e0;
OSUM = 0.0e0;
@ -566,7 +566,7 @@ sit(void)
C
*/
AGAMMA = 3*sit_A0; /* Grenthe p 379 */
A = AGAMMA / log(10.0);
A = AGAMMA / LOG_10;
/*
* F is now for log10 gamma
*/
@ -652,7 +652,7 @@ sit(void)
C
*/
/*COSMOT = 1.0e0 + 2.0e0 * OSMOT / OSUM;*/
COSMOT = 1.0e0 + OSMOT*log(10.0) / OSUM;
COSMOT = 1.0e0 + OSMOT*LOG_10 / OSUM;
/*
C
C CALCULATE THE ACTIVITY OF WATER
@ -1013,7 +1013,7 @@ Restart:
base[i] = residual[i];
}
d = 0.0001;
d1 = d * log(10.0);
d1 = d * LOG_10;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
{
@ -1096,7 +1096,7 @@ Restart:
residuals();
for (j = 0; j < count_unknowns; j++)
{
array[j * (count_unknowns + 1) + i] =
my_array[j * (count_unknowns + 1) + i] =
-(residual[j] - base[j]) / d2;
}
switch (x[i]->type)
@ -1115,9 +1115,9 @@ Restart:
break;
case MH:
s_eminus->la -= d;
if (array[i * (count_unknowns + 1) + i] == 0)
if (my_array[i * (count_unknowns + 1) + i] == 0)
{
array[i * (count_unknowns + 1) + i] =
my_array[i * (count_unknowns + 1) + i] =
exp(s_h2->lm * LOG_10) * 2;
}
break;
@ -1227,6 +1227,7 @@ model_sit(void)
PhreeqcIWait(this);
#endif
iterations++;
overall_iterations++;
if (iterations > itmax - 1 && debug_model == FALSE
&& pr.logfile == TRUE)
{

27
src/phreeqcpp/tidy.cpp
View File

@ -20,6 +20,7 @@ tidy_model(void)
/*
* Determine if any new elements, species, phases have been read
*/
overall_iterations = 0;
state = INITIALIZE;
new_model = FALSE;
new_pp_assemblage = FALSE;
@ -4005,16 +4006,16 @@ tidy_min_surface(void)
free_check_null(temp_formula);
continue;
}
if (strcmp(elt_ptr->master->s->name, temp_formula) != 0)
{
error_string = sformatf("Suggest using master species formula in SURFACE \n\t for surface related to equilibrium_phase: %s.",
elt_ptr->master->s->name);
warning_msg(error_string);
}
if (elt_ptr->master->s->z != 0.0)
//if (strcmp(elt_ptr->master->s->name, temp_formula) != 0)
//{
// error_string = sformatf("Suggest using master species formula in SURFACE \n\t for surface related to equilibrium_phase: %s.",
// elt_ptr->master->s->name);
// warning_msg(error_string);
//}
if (elt_ptr->master->s->z != 0.0 && surface_ptr->Get_dl_type() != cxxSurface::DONNAN_DL)
{
error_string = sformatf(
"Suggest master species of surface, %s, be uncharged for surface related to equilibrium_phase.",
"Use the -donnan option when coupling surface %s to an equilibrium_phase, \n\t and note to give the equilibrium_phase the surface charge.",
elt_ptr->master->s->name);
warning_msg(error_string);
}
@ -4032,6 +4033,16 @@ tidy_min_surface(void)
Further, if you precipitate Ca-Mont, make SurfCa, desorb
all the Ca, then dissolve the "Ca-Mont", you must remove SurfCa, or you
will end up with Ca in solution. H and O are excluded */
/* Example that makes montmorillonite a cation exchanger:
PHASES
Summ_Montmorillonite; Al2.33Si3.67O10(OH)2-0.33 + 12 H2O = 2.33 Al(OH)4- + 3.67 H4SiO4 + 2 H+; -log_k -44.4
SURFACE_MASTER_SPECIES; Summ Summ-; SURFACE_SPECIES; Summ- = Summ-
SOLUTION 1; Na 1e1; Cl 1e1; pH 7 charge; C(4) 1 CO2(g) -2
EQUILIBRIUM_PHASES 1; Ca-Montmorillonite 0 1e-3
Summ_Montmorillonite 0 0
SURFACE 1; Summ Summ_Montmorillonite 0.33 3.11e5; -donnan; -equil 1
END
*/
for (int jj = 0; jj < count_elts; jj++)
{
if (elt_list[jj].elt->primary->s->type != SURF

71
src/phreeqcpp/transport.cpp
View File

@ -431,7 +431,7 @@ transport(void)
/*
* Now transport
*/
sprintf(token, "\nCalculating transport: %d cells, %d shifts, %d mixruns...\n\n",
sprintf(token, "\nCalculating transport: %d (mobile) cells, %d shifts, %d mixruns...\n\n",
count_cells, count_shifts - transport_start + 1, nmix);
screen_msg(token);
max_iter = 0;
@ -514,8 +514,8 @@ transport(void)
{
if (!dV_dcell && (i == 0 || i == count_cells + 1))
continue;
if (iterations > max_iter)
max_iter = iterations;
if (overall_iterations > max_iter)
max_iter = overall_iterations;
cell_no = i;
mixrun = j;
if (multi_Dflag)
@ -708,8 +708,8 @@ transport(void)
run_reactions(i, kin_time, NOMIX, step_fraction);
if (multi_Dflag == TRUE)
fill_spec(i);
if (iterations > max_iter)
max_iter = iterations;
if (overall_iterations > max_iter)
max_iter = overall_iterations;
if (nmix == 0 && stag_data->count_stag == 0)
print_punch(i, true);
if (i == first_c && count_cells > 1)
@ -811,8 +811,8 @@ transport(void)
{
if (!dV_dcell && (i == 0 || i == count_cells + 1))
continue;
if (iterations > max_iter)
max_iter = iterations;
if (overall_iterations > max_iter)
max_iter = overall_iterations;
cell_no = i;
mixrun = j;
if (multi_Dflag)
@ -973,8 +973,8 @@ void Phreeqc::
print_punch(int i, boolean active)
/* ---------------------------------------------------------------------- */
{
if ((!(cell_data[i].punch && (transport_step % punch_modulus == 0)) &&
!(cell_data[i].print && (transport_step % print_modulus == 0))) ||
if (!(cell_data[i].punch && (transport_step % punch_modulus == 0)) &&
!(cell_data[i].print && (transport_step % print_modulus == 0)) ||
(bcon_first == 2 && i == 0) ||
(bcon_last == 2 && i == count_cells + 1))
return;
@ -2510,10 +2510,6 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
else if (icell == count_cells)
ct[icell].visc2 = ct[icell].visc1;
}
//LDBLE d_damper = Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_mu() /
// Utilities::Rxn_find(Rxn_solution_map, icell)->Get_mu();
//d_damper = pow(d_damper, 0.3);
//if (d_damper > 1) ct[icell].visc1 *= d_damper; else ct[icell].visc2 *= d_damper;
/* in each cell: DL surface = mass_water_DL / (cell_length)
free pore surface = mass_water_free / (cell_length)
@ -2757,7 +2753,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
if (ct[icell].v_m[k].z)
ct[icell].v_m[k].zc = ct[icell].v_m[k].z * c1;
if (dV_dcell && ct[icell].v_m[k].z)
if (dV_dcell && ct[icell].v_m[k].z && !fix_current)
{
// compare diffusive and electromotive forces
dum = ct[icell].v_m[k].grad;
@ -2817,6 +2813,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
b_j = A2 * (f_free_j + g_j / ct[icell].visc2);
if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
else if (icell == all_cells - 1 && stagnant)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf *= 2 for this 'reservoir' cell in the input */
else
{
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
@ -2828,15 +2826,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
dum2 *= sol_D[jcell].viscos_f;
b_j *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == 0 && !stagnant)
{
//if (d_damper < 1 && g_i == 0)
// ct[icell].v_m[k].b_ij = A2 * sol_D[icell].spec[i].Dwt * (f_free_j + g_j * d_damper / ct[icell].visc2);
//else
ct[icell].v_m[k].b_ij = b_j;
}
else
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
else if (icell == 3 && stagnant && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for stagnant cell 3 in the input */
}
if (ct[icell].v_m[k].z)
@ -2874,7 +2868,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
if (ct[icell].v_m[k].z)
ct[icell].v_m[k].zc = ct[icell].v_m[k].z * c2;
if (dV_dcell && ct[icell].v_m[k].z)
if (dV_dcell && ct[icell].v_m[k].z && !fix_current)
{
// compare diffuse and electromotive forces
dum = ct[icell].v_m[k].grad;
@ -2933,6 +2927,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
if (icell == 0 && !stagnant)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == 3 && stagnant && g_j && !g_i)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for 'reservoir' cell 3 in the input */
else
{
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
@ -2944,10 +2940,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
dum2 *= sol_D[icell].viscos_f;
b_i *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
else
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
else if (jcell == all_cells - 1 && stagnant && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf * 2 for this 'reservoir' cell in the input */
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
@ -2989,7 +2986,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
if (ct[icell].v_m[k].z)
ct[icell].v_m[k].zc = ct[icell].v_m[k].z * ct[icell].v_m[k].c;
if (dV_dcell && ct[icell].v_m[k].z)
if (dV_dcell && ct[icell].v_m[k].z && !fix_current)
{
// compare diffuse and electromotive forces
dum = ct[icell].v_m[k].grad;
@ -3035,13 +3032,23 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
if (icell == 0 && !stagnant)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
// but for boundary cells...
if (stagnant)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2;
}
else
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
{
if (icell == 0)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == count_cells)
ct[icell].v_m[k].b_ij = b_i;
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
@ -3067,6 +3074,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
if (dV_dcell)
{
if (transport_step >= 100)
icell = icell;
current_cells[icell].ele = current_cells[icell].dif = 0;
dum = dV_dcell * F_Re3 / tk_x2;
for (i = 0; i < ct[icell].J_ij_count_spec; i++)

6
src/phreeqcpp/utilities.cpp
View File

@ -172,8 +172,10 @@ calc_rho_0(LDBLE tc, LDBLE pa)
if (need_temp_msg < 1)
{
std::ostringstream w_msg;
w_msg << "Fitting range for density of pure water is 0-300 C.\n";
w_msg << "Using temperature of 350 C for density and dielectric calculation.";
w_msg << "Fitting range for dielectric constant of pure water is 0-350 C.\n";
w_msg << "Fitting range for density along the saturation pressure line is 0-374 C,\n";
w_msg << " for higher pressures up to 1000 atm 0-300 C.\n";
w_msg << "Using temperature of 350 C for dielectric and density calculation.";
warning_msg(w_msg.str().c_str());
need_temp_msg++;
}