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 93ab5c9..82515f7

Browse Source 
82515f7 Merge remote-tracking branch 'origin/master'
1cc4fa0 fix phi for water at high TP. Tony May, 2023
275b340 Tony's changes May 5 and 7.
0fb2509 cleaned up compiler warnings. removed math.h. cleaned up NAN

git-subtree-dir: phreeqcpp
git-subtree-split: 82515f7fecada9be373fd832dd0c5cdf53c79da7
This commit is contained in:
Darth Vader 2023-05-26 04:27:58 +00:00
parent dc26aa2a8f
commit 050663209b
19 changed files with 78 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
ChartHandler.cpp
View File

@ -170,18 +170,14 @@ ChartHandler::End_timer()
io->error_flush();
}
}
size_t i(0), i2(0);
for ( ; it != chart_map.end(); it++)
{
i = 0;
it->second->Rate_free();
if (it->second->Get_form_started())
{
#if defined(__cplusplus_cli)
while (0 != System::Threading::Interlocked::CompareExchange(it->second->usingResource, 6, 0))
{
//if (i > max_tries) break;
i++;
System::Threading::Thread::Sleep(60);
}
#endif
@ -191,12 +187,8 @@ ChartHandler::End_timer()
int n = System::Threading::Interlocked::Exchange(it->second->usingResource, 0);
assert(n == 6);
#endif
i2 = 0;
while (it->second->Get_done() != true)
{
//if (i2 > max_tries) break;
i2++;
#if defined(__cplusplus_cli)
System::Threading::Thread::Sleep(60);
#endif
@ -225,7 +217,6 @@ ChartHandler::dump(std::ostream & oss, unsigned int indent)
std::map<int, ChartObject *>::iterator it = this->chart_map.begin();
for ( ; it != chart_map.end(); it++)
{
size_t i = 0;
it->second->dump(oss, indent);
}
return true;

1
ChartObject.cpp
View File

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

11
PBasic.cpp
View File

@ -1558,6 +1558,9 @@ listtokens(FILE * f, tokenrec * l_buf)
case tokt_sc:
output_msg("T_SC");
break;
case tokf_visc:
output_msg("F_VISC");
break;
case toktc:
output_msg("TC");
break;
@ -3904,6 +3907,13 @@ factor(struct LOC_exec * LINK)
}
break;
case tokf_visc:
{
const char* str = stringfactor(STR1, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->calc_f_visc(str);
}
break;
case toktc:
{
n.UU.val = PhreeqcPtr->tc_x;
@ -7518,6 +7528,7 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("surf", PBasic::toksurf),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sys", PBasic::toksys),
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("f_visc", PBasic::tokf_visc),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("tc", PBasic::toktc),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("time", PBasic::toktime),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("title", PBasic::toktitle),

2
PBasic.h
View File

@ -7,7 +7,6 @@
#include <stdio.h>
#include <limits.h>
#include <ctype.h>
//#include <math.h>
#include <cmath>
#include <setjmp.h>
#include "phrqtype.h"
@ -338,6 +337,7 @@ public:
toktotmol,
toktotmoles,
toktrim,
tokf_visc,
tokviscos,
tokviscos_0,
tokvm,

1
Phreeqc.h
View File

@ -996,6 +996,7 @@ public:
int reformat_surf(const char* comp_name, LDBLE fraction, const char* new_comp_name,
LDBLE new_Dw, int cell);
LDBLE viscosity(void);
LDBLE calc_f_visc(const char *name);
LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);

5
SolutionIsotope.cxx
View File

@ -67,12 +67,7 @@ cxxSolutionIsotope::dump_xml(std::ostream & s_oss, unsigned int indent) const
s_oss << indent1;
s_oss << "iso_ratio=\"" << this->ratio << "\"" << "\n";
#ifdef NPP
if (!isnan(this->ratio_uncertainty))
#else
//if (this->ratio_uncertainty != NAN)
if (!std::isnan(this->ratio_uncertainty))
#endif
{
s_oss << indent1;
s_oss << "iso_ratio_uncertainty=\"" << this->

64
basicsubs.cpp
View File

@ -1252,6 +1252,24 @@ calc_t_sc(const char *name)
return (-999.99);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_f_visc(const char *name)
/* ---------------------------------------------------------------------- */
{
char token[MAX_LENGTH];
class species *s_ptr;
if (print_viscosity)
{
strcpy(token, name);
s_ptr = s_search(token);
if (s_ptr != NULL)
return s_ptr->dw_t_visc;
}
return 0;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
equi_phase(const char *phase_name)
@ -1599,12 +1617,7 @@ get_calculate_value(const char *name)
calculate_value_ptr->name);
error_msg(error_string, STOP);
}
#ifdef NPP
if (isnan(rate_moles))
#else
//if (rate_moles == NAN)
if(std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
calculate_value_ptr->name);
@ -1762,18 +1775,17 @@ LDBLE Phreeqc::
pr_pressure(const char* phase_name)
/* ---------------------------------------------------------------------- */
{
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);
}
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
if (gas_phase_ptr != NULL)
{
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);
}
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]);
@ -1792,6 +1804,10 @@ pr_pressure(const char* phase_name)
}
}
}
else if (phase_ptr->in != FALSE && phase_ptr->pr_in)
{
return phase_ptr->pr_p;
}
return(0.0);
}
/* ---------------------------------------------------------------------- */
@ -1807,35 +1823,35 @@ LDBLE Phreeqc::
pr_phi(const char *phase_name)
/* ---------------------------------------------------------------------- */
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
if (gas_phase_ptr != NULL)
{
int l;
class phase* phase_ptr = phase_bsearch(phase_name, &l, FALSE);
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);
}
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
if (gas_phase_ptr != NULL)
{
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);
if (phase_ptr == phase_ptr_gas)
{
if (gas_phase_ptr->Get_pr_in())
{
return phase_ptr->pr_phi;
}
return phase_ptr->pr_phi;
else
{
return gas_comp_ptr->Get_phi();
}
}
}
}
else if (phase_ptr->in != FALSE && phase_ptr->pr_in)
{
return phase_ptr->pr_phi;
}
return (1.0);
}
/* ---------------------------------------------------------------------- */

1
cl1.cpp
View File

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

1
common/Utils.cxx
View File

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

1
gases.cpp
View File

@ -607,6 +607,7 @@ calc_PR(void)
phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi));
}
else
phi = -3.0; // fugacity coefficient = 0.05

2
global_structures.h
View File

@ -715,6 +715,7 @@ public:
dw_a2 = 0;
dw_a_visc = 0; // viscosity correction of SC
dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC
dw_t_visc = 0; // contribution to viscosity
dw_corr = 0; // dw corrected for TK and mu
erm_ddl = 0; // enrichment factor in DDL
equiv = 0; // equivalents in exchange species
@ -776,6 +777,7 @@ public:
LDBLE dw_a2;
LDBLE dw_a_visc;
LDBLE dw_t_SC;
LDBLE dw_t_visc;
LDBLE dw_corr;
LDBLE erm_ddl;
LDBLE equiv;

24
inverse.cpp
View File

@ -3582,40 +3582,21 @@ check_isotopes(class inverse *inv_ptr)
i = ii;
/* use inverse-defined uncertainties first */
#ifdef NPP
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()
&& !std::isnan(inv_ptr->i_u[i].uncertainties[j]))
#endif
{
kit->second.Set_x_ratio_uncertainty(inv_ptr->i_u[i].uncertainties[j]);
/* use solution-defined uncertainties second */
}
#ifdef NPP
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
&& !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]);
/* use solution-defined uncertainties second */
}
#ifdef NPP
else if (!isnan(kit->second.Get_ratio_uncertainty()))
#else
//else if (kit->second.Get_ratio_uncertainty() != NAN)
else if (!std::isnan(kit->second.Get_ratio_uncertainty()))
#endif
{
kit->second.Set_x_ratio_uncertainty(
kit->second.Get_ratio_uncertainty());
@ -3643,12 +3624,7 @@ 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 (std::isnan(kit->second.Get_x_ratio_uncertainty()))
#endif
{
error_string = sformatf(
"In solution %d, isotope ratio uncertainty is needed for element: %g%s.",

15
isotopes.cpp
View File

@ -900,12 +900,7 @@ punch_calculate_values(void)
calculate_value_ptr->name);
error_msg(error_string, STOP);
}
#ifdef NPP
if (isnan(rate_moles))
#else
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
calculate_value_ptr->name);
@ -1134,12 +1129,7 @@ calculate_values(void)
calculate_value_ptr->name);
error_msg(error_string, STOP);
}
#ifdef NPP
if (isnan(rate_moles))
#else
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
calculate_value_ptr->name);
@ -1202,12 +1192,7 @@ calculate_values(void)
calculate_value_ptr->name);
error_msg(error_string, STOP);
}
#ifdef NPP
if (isnan(rate_moles))
#else
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Calculated value not SAVEed for %s.",
calculate_value_ptr->name);

5
kinetics.cpp
View File

@ -115,12 +115,7 @@ calc_kinetic_reaction(cxxKinetics *kinetics_ptr, LDBLE time_step)
kinetics_comp_ptr->Get_rate_name().c_str());
error_msg(error_string, STOP);
}
#ifdef NPP
if (isnan(rate_moles))
#else
//if (rate_moles == NAN)
if (std::isnan(rate_moles))
#endif
{
error_string = sformatf( "Moles of reaction not SAVEed for %s.",
kinetics_comp_ptr->Get_rate_name().c_str());

1
prep.cpp
View File

@ -3997,6 +3997,7 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
phi = (phi > 4.44 ? 4.44 : (phi < -4.6 ? -4.6 : phi));
//if (phi > 4.44)
// phi = 4.44;
}

8
read.cpp
View File

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

1
sundialsmath.cpp
View File

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

20
tidy.cpp
View File

@ -979,12 +979,7 @@ tidy_gas_phase(void)
error_msg(error_string, CONTINUE);
}
/* calculate moles */
#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 (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#endif
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
if (!PR)
@ -1013,12 +1008,7 @@ tidy_gas_phase(void)
*/
if (!gas_phase_ptr->Get_solution_equilibria())
{
#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 (!std::isnan(gas_phase_ptr->Get_gas_comps()[j].Get_p_read()))
#endif
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
if (!PR)
@ -1689,12 +1679,7 @@ tidy_ss_assemblage(void)
phase_ptr->moles_x = 0;
phase_ptr->fraction_x = 0;
}
#ifdef NPP
if (isnan(comp_ptr->Get_moles()))
#else
//if (comp_ptr->Get_moles() == NAN)
if (std::isnan(comp_ptr->Get_moles()))
#endif
{
input_error++;
error_string = sformatf(
@ -3021,12 +3006,7 @@ tidy_isotopes(void)
temp_iso.Set_total(0);
temp_iso.Set_ratio(master[k]->isotope_ratio);
temp_iso.Set_ratio_uncertainty(master[k]->isotope_ratio_uncertainty);
#ifdef NPP
if (!isnan(master[k]->isotope_ratio_uncertainty))
#else
//if (master[k]->isotope_ratio_uncertainty != NAN)
if (!std::isnan(master[k]->isotope_ratio_uncertainty))
#endif
{
temp_iso.Set_ratio_uncertainty_defined(true);
}

31
transport.cpp
View File

@ -5963,10 +5963,10 @@ viscosity(void)
We use the harmonic mean of the Dw's, and the arithmetic mean of the z's,
both weighted by the equivalent concentration.
*/
LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, ta;
LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, t3, fan = 1;
LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, tc;
m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = ta = 0;
m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = 0;
tc = (tc_x > 200) ? 200 : tc_x;
@ -5978,6 +5978,7 @@ viscosity(void)
continue;
if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3])
{
s_x[i]->dw_t_visc = 0;
t1 = s_x[i]->moles / mass_water_aq_x;
l_z = fabs(s_x[i]->z);
if (l_z)
@ -5993,8 +5994,9 @@ viscosity(void)
s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0);
}
// find B * m and D * m * mu^d3
Bc += (s_x[i]->Jones_Dole[0] +
s_x[i]->dw_t_visc = (s_x[i]->Jones_Dole[0] +
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
Bc += s_x[i]->dw_t_visc;
// 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];
@ -6002,8 +6004,10 @@ viscosity(void)
t2 = -0.8 / s_x[i]->Jones_Dole[5];
else
t2 = -1;
Dc += (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) *
t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) *
t1 * (pow(mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2);
s_x[i]->dw_t_visc += t3;
Dc += t3;
//output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc ));
}
// parms for A...
@ -6023,13 +6027,11 @@ viscosity(void)
{
V_Cl = s_x[i]->logk[vm_tc];
V_an += V_Cl * s_x[i]->moles;
ta += s_x[i]->moles;
m_an += s_x[i]->moles;
}
else// if (s_x[i]->Jones_Dole[6])
{
V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * s_x[i]->moles;
ta += s_x[i]->moles;
m_an += s_x[i]->moles;
}
if (Dw)
@ -6064,23 +6066,24 @@ viscosity(void)
A = 0;
viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / mass_water_aq_x);
if (m_an)
{
V_an /= m_an;
ta /= m_an;
}
if (!V_Cl)
V_Cl = calc_vm_Cl();
if (V_an && V_Cl && ta)
viscos += (viscos_0 * (2 - ta * V_an / V_Cl) * (Bc + Dc));
else
viscos += (viscos_0 * (Bc + Dc));
if (V_an)
fan = 2 - V_an / V_Cl;
//else
// fan = 1;
viscos += viscos_0 * fan * (Bc + Dc);
if (viscos < 0)
{
viscos = viscos_0;
warning_msg("viscosity < 0, reset to viscosity of pure water");
}
for (i = 0; i < (int)this->s_x.size(); i++)
{
s_x[i]->dw_t_visc /= (Bc + Dc);
}
return viscos;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::