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 '45befa710880429a2a04cbc0294e83e0f2ce66b9'

Browse Source
This commit is contained in:
Darth Vader 2024-08-29 19:02:04 +00:00
commit 7e60110c31
5 changed files with 174 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
phreeqcpp/PBasic.cpp
View File

@ -4732,6 +4732,8 @@ factor(struct LOC_exec * LINK)
case tokviscos: case tokviscos:
{ {
if (PhreeqcPtr->print_viscosity)
PhreeqcPtr->viscosity(nullptr);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->viscos; n.UU.val = (parse_all) ? 1 : PhreeqcPtr->viscos;
} }
break; break;

37
phreeqcpp/basicsubs.cpp
View File

@ -205,8 +205,11 @@ diff_c(const char *species_name)
Dw *= viscos_0_25 / viscos_0; Dw *= viscos_0_25 / viscos_0;
} }
if (s_ptr->dw_a_v_dif) if (s_ptr->dw_a_v_dif && print_viscosity)
{
viscosity(nullptr);
Dw *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif); Dw *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif);
}
return Dw; return Dw;
} }
@ -235,8 +238,11 @@ setdiff_c(const char *species_name, double d, double d_v_d)
Dw *= viscos_0_25 / viscos_0; Dw *= viscos_0_25 / viscos_0;
} }
if (d_v_d) if (d_v_d && print_viscosity)
Dw *= pow(viscos_0 / viscos, d_v_d); {
viscosity(nullptr);
Dw *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif);
}
return Dw; return Dw;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -269,6 +275,8 @@ calc_SC(void)
// } // }
//} //}
av = 0; av = 0;
if (print_viscosity)
viscosity(nullptr);
if (!Falk) if (!Falk)
{ {
for (i = 0; i < (int)this->s_x.size(); i++) for (i = 0; i < (int)this->s_x.size(); i++)
@ -1171,6 +1179,29 @@ diff_layer_total(const char* total_name, const char* surface_name)
return (0); return (0);
} }
} }
else if (strcmp_nocase("viscos_ddl", total_name) == 0)
{
if (dl_type_x != cxxSurface::NO_DL)
{
cxxSurfaceCharge* charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
if (charge_ptr->Get_mass_water() > 0)
{
cxxSurface * surf_ptr = use.Get_surface_ptr();
if (surf_ptr->Get_calc_viscosity())
{
viscosity(surf_ptr);
viscosity(nullptr);
return charge_ptr->Get_DDL_viscosity();
}
else
return charge_ptr->Get_DDL_viscosity() * viscos;
}
}
else
{
return (0);
}
}
/* /*
* find total moles of each element in diffuse layer... * find total moles of each element in diffuse layer...
*/ */

85
phreeqcpp/class_main.cpp
View File

@ -1,6 +1,8 @@
#ifndef NPP
#ifdef DOS #ifdef DOS
#include <windows.h> #include <windows.h>
#endif #endif
#endif
#include "Phreeqc.h" #include "Phreeqc.h"
#include "NameDouble.h" #include "NameDouble.h"
@ -122,12 +124,18 @@ main_method(int argc, char *argv[])
{ {
return errors; return errors;
} }
#ifdef NPP
#ifdef DOS
write_banner();
#endif
#else
#ifndef NO_UTF8_ENCODING #ifndef NO_UTF8_ENCODING
#ifdef DOS #ifdef DOS
SetConsoleOutputCP(CP_UTF8); SetConsoleOutputCP(CP_UTF8);
#endif #endif
write_banner(); write_banner();
#endif #endif
#endif
/* /*
* Initialize arrays * Initialize arrays
@ -214,11 +222,17 @@ main_method(int argc, char *argv[])
{ {
return errors; return errors;
} }
#ifdef NPP
#ifdef DOS
write_banner();
#endif
#else
#ifndef NO_UTF8_ENCODING #ifndef NO_UTF8_ENCODING
#ifdef DOS #ifdef DOS
SetConsoleOutputCP(CP_UTF8); SetConsoleOutputCP(CP_UTF8);
#endif #endif
write_banner(); write_banner();
#endif
#endif #endif
/* /*
@ -281,6 +295,64 @@ main_method(int argc, char *argv[])
return 0; return 0;
} }
#endif //TEST_COPY #endif //TEST_COPY
#ifdef NPP
/* ---------------------------------------------------------------------- */
int Phreeqc::
write_banner(void)
/* ---------------------------------------------------------------------- */
{
#ifdef TESTING
return OK;
#endif
#ifndef NO_UTF8_ENCODING
char buffer[80];
int len, indent;
screen_msg(
" ÛßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÛ\n");
screen_msg(
" º º\n");
/* version */
#ifdef NPP
len = sprintf(buffer, "* PHREEQC-%s *", "3.8.0");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
indent = (49 - len) / 2;
screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer,
49 - indent - len, ' '));
screen_msg(
" º º\n");
screen_msg(
" º A hydrogeochemical transport model º\n");
screen_msg(
" º º\n");
screen_msg(
" º by º\n");
screen_msg(
" º D.L. Parkhurst and C.A.J. Appelo º\n");
screen_msg(
" º º\n");
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "August 27, 2024, with bug-fixes and new items");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
indent = (49 - len) / 2;
screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer,
49 - indent - len, ' '));
screen_msg(
" ÛÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÛ\n\n");
#endif
return 0;
}
#else
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int Phreeqc:: int Phreeqc::
write_banner(void) write_banner(void)
@ -298,11 +370,7 @@ write_banner(void)
" ║ ║\n"); " ║ ║\n");
/* version */ /* version */
#ifdef NPP
len = sprintf(buffer, "* PHREEQC-%s *", "3.7.3");
#else
len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "@VERSION@"); len = snprintf(buffer, sizeof(buffer), "* PHREEQC-%s *", "@VERSION@");
#endif
indent = (44 - len) / 2; indent = (44 - len) / 2;
screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer, screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer,
44 - indent - len, ' ')); 44 - indent - len, ' '));
@ -322,11 +390,7 @@ write_banner(void)
/* date */ /* date */
#ifdef NPP
len = sprintf(buffer, "%s", "March 14, 2024, with bug-fixes and new items");
#else
len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@"); len = snprintf(buffer, sizeof(buffer), "%s", "@VER_DATE@");
#endif
indent = (44 - len) / 2; indent = (44 - len) / 2;
screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer, screen_msg(sformatf("%14c║%*c%s%*c║\n", ' ', indent, ' ', buffer,
44 - indent - len, ' ')); 44 - indent - len, ' '));
@ -336,6 +400,7 @@ write_banner(void)
#endif #endif
return 0; return 0;
} }
#endif
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int Phreeqc:: int Phreeqc::
@ -507,7 +572,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
output_msg(sformatf(" Input file: %s\n", in_file.c_str())); output_msg(sformatf(" Input file: %s\n", in_file.c_str()));
output_msg(sformatf(" Output file: %s\n", out_file.c_str())); output_msg(sformatf(" Output file: %s\n", out_file.c_str()));
#ifdef NPP #ifdef NPP
output_msg(sformatf("Using PHREEQC: version 3.7.3, compiled March 14, 2024, with bug-fixes and new items\n")); output_msg(sformatf("Using PHREEQC: version 3.8.2, compiled August 27, 2024, with bug-fixes and new items\n"));
#endif #endif
output_msg(sformatf("Database file: %s\n\n", token.c_str())); output_msg(sformatf("Database file: %s\n\n", token.c_str()));
#ifdef NPP #ifdef NPP
@ -516,7 +581,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
/* /*
* local cleanup * local cleanup
*/ */
line = (char *) free_check_null(line); line = (char *) free_check_null(line);
line_save = (char *) free_check_null(line_save); line_save = (char *) free_check_null(line_save);

26
phreeqcpp/print.cpp
View File

@ -1,3 +1,4 @@
//Note to encode in ANSI with NP++
#include "Utils.h" #include "Utils.h"
#include "Phreeqc.h" #include "Phreeqc.h"
#include "phqalloc.h" #include "phqalloc.h"
@ -270,14 +271,23 @@ print_diffuse_layer(cxxSurfaceCharge *charge_ptr)
output_msg(sformatf( output_msg(sformatf(
"\tWater in diffuse layer: %8.3e kg, %4.1f%% of total DDL-water.\n", "\tWater in diffuse layer: %8.3e kg, %4.1f%% of total DDL-water.\n",
(double) charge_ptr->Get_mass_water(), (double) d)); (double) charge_ptr->Get_mass_water(), (double) d));
if (charge_ptr->Get_DDL_viscosity()) if (print_viscosity && d > 0)
{ {
if (d == 100) cxxSurface * surf_ptr = use.Get_surface_ptr();
output_msg(sformatf( if (surf_ptr->Get_calc_viscosity())
"\t\t viscosity: %7.5f mPa s.\n", (double)charge_ptr->Get_DDL_viscosity())); {
viscosity(surf_ptr);
viscosity(nullptr);
if (d == 100)
output_msg(sformatf(
"\t\t calculated viscosity: %7.5f mPa s.\n", (double)charge_ptr->Get_DDL_viscosity()));
else
output_msg(sformatf(
"\t\t calculated viscosity: %7.5f mPa s for this DDL water. (%7.5f mPa s for total DDL-water.)\n", (double)charge_ptr->Get_DDL_viscosity(), (double)use.Get_surface_ptr()->Get_DDL_viscosity()));
}
else else
output_msg(sformatf( output_msg(sformatf(
"\t\t viscosity: %7.5f mPa s for this DDL water. (%7.5f mPa s for total DDL-water.)\n", (double)charge_ptr->Get_DDL_viscosity(), (double)use.Get_surface_ptr()->Get_DDL_viscosity())); "\t\t viscosity: %7.5f mPa s for DDL water.\n", (double)charge_ptr->Get_DDL_viscosity() * viscos));
} }
if (use.Get_surface_ptr()->Get_debye_lengths() > 0 && d > 0) if (use.Get_surface_ptr()->Get_debye_lengths() > 0 && d > 0)
@ -2256,8 +2266,8 @@ print_totals(void)
//#ifdef NPP //#ifdef NPP
if (print_viscosity) if (print_viscosity)
{ {
output_msg(sformatf("%45s%9.5f", "Viscosity (mPa s) = ", viscosity(nullptr);
(double) viscos)); output_msg(sformatf("%45s%9.5f", "Viscosity (mPa s) = ", (double) viscos));
if (tc_x > 200 && !pure_water) if (tc_x > 200 && !pure_water)
{ {
#ifdef NO_UTF8_ENCODING #ifdef NO_UTF8_ENCODING

71
phreeqcpp/transport.cpp
View File

@ -962,6 +962,7 @@ transport(void)
{ {
snprintf(token, sizeof(token), snprintf(token, sizeof(token),
"\nFor balancing negative concentrations in MCD, added in total to the system:"); "\nFor balancing negative concentrations in MCD, added in total to the system:");
count_warnings = pr.warnings - 1;
warning_msg(token); warning_msg(token);
for (i = 0; i < count_moles_added; i++) for (i = 0; i < count_moles_added; i++)
{ {
@ -970,6 +971,7 @@ transport(void)
snprintf(token, sizeof(token), snprintf(token, sizeof(token),
"\t %.4e moles %s.", "\t %.4e moles %s.",
(double)moles_added[i].moles, moles_added[i].name); (double)moles_added[i].moles, moles_added[i].name);
count_warnings = pr.warnings - 1;
warning_msg(token); warning_msg(token);
} }
} }
@ -3937,15 +3939,20 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{ {
if (s_ptr1->Get_dl_type() != cxxSurface::NO_DL) if (s_ptr1->Get_dl_type() != cxxSurface::NO_DL)
{ {
if (s_ptr1->Get_calc_viscosity())
{
viscosity(s_ptr1);
ct[icell].visc1 = s_ptr1->Get_DDL_viscosity();
}
else
ct[icell].visc1 = s_ptr1->Get_DDL_viscosity() * Utilities::Rxn_find(Rxn_solution_map, icell)->Get_viscosity();
s_charge_p.assign(s_ptr1->Get_surface_charges().begin(), s_ptr1->Get_surface_charges().end()); s_charge_p.assign(s_ptr1->Get_surface_charges().begin(), s_ptr1->Get_surface_charges().end());
s_com_p.assign(s_ptr1->Get_surface_comps().begin(), s_ptr1->Get_surface_comps().end()); s_com_p.assign(s_ptr1->Get_surface_comps().begin(), s_ptr1->Get_surface_comps().end());
if (s_ptr1->Get_only_counter_ions()) if (s_ptr1->Get_only_counter_ions())
only_counter = TRUE; only_counter = TRUE;
ct[icell].visc1 = s_ptr1->Get_DDL_viscosity();
if (s_ptr1->Get_calc_viscosity())
ct[icell].visc1 /= Utilities::Rxn_find(Rxn_solution_map, icell)->Get_viscosity();
/* find the immobile surface charges with DL... */ /* find the immobile surface charges with DL... */
for (i = 0; i < (int)s_charge_p.size(); i++) for (i = 0; i < (int)s_charge_p.size(); i++)
{ {
@ -3966,16 +3973,20 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{ {
if (s_ptr2->Get_dl_type() != cxxSurface::NO_DL) if (s_ptr2->Get_dl_type() != cxxSurface::NO_DL)
{ {
if (s_ptr2->Get_calc_viscosity())
{
viscosity(s_ptr2);
ct[icell].visc2 = s_ptr2->Get_DDL_viscosity();
}
else
ct[icell].visc2 = s_ptr2->Get_DDL_viscosity() * Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_viscosity();
s_charge_p.assign(s_ptr2->Get_surface_charges().begin(), s_ptr2->Get_surface_charges().end()); s_charge_p.assign(s_ptr2->Get_surface_charges().begin(), s_ptr2->Get_surface_charges().end());
s_com_p.assign(s_ptr2->Get_surface_comps().begin(), s_ptr2->Get_surface_comps().end()); s_com_p.assign(s_ptr2->Get_surface_comps().begin(), s_ptr2->Get_surface_comps().end());
if (s_ptr2->Get_only_counter_ions()) if (s_ptr2->Get_only_counter_ions())
only_counter = TRUE; only_counter = TRUE;
ct[icell].visc2 = s_ptr2->Get_DDL_viscosity();
if (s_ptr2->Get_calc_viscosity())
ct[icell].visc2 /= Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_viscosity();
for (i = 0; i < (int)s_charge_p.size(); i++) for (i = 0; i < (int)s_charge_p.size(); i++)
{ {
for (i1 = 0; i1 < (int)s_com_p.size(); i1++) for (i1 = 0; i1 < (int)s_com_p.size(); i1++)
@ -3990,6 +4001,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
} }
} }
} }
viscosity(nullptr);
if (!stagnant) if (!stagnant)
{ {
if (icell == 0) if (icell == 0)
@ -5975,8 +5987,11 @@ viscosity(cxxSurface *surf_ptr)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
if (surf_ptr && !surf_ptr->Get_calc_viscosity()) if (surf_ptr && !surf_ptr->Get_calc_viscosity())
{
for (int i = 0; i < (int)surf_ptr->Get_surface_charges().size(); i++)
surf_ptr->Get_surface_charges()[i].Set_DDL_viscosity(surf_ptr->Get_DDL_viscosity());
return surf_ptr->Get_DDL_viscosity(); return surf_ptr->Get_DDL_viscosity();
}
/* from Atkins, 1994. Physical Chemistry, 5th ed. */ /* from Atkins, 1994. Physical Chemistry, 5th ed. */
//viscos = //viscos =
@ -6120,7 +6135,7 @@ viscosity(cxxSurface *surf_ptr)
eq_plus += l_moles; eq_plus += l_moles;
l_mu_x += l_moles * s_x[i]->z; l_mu_x += l_moles * s_x[i]->z;
} }
l_mu_x += fabs(eq_plus + eq_min); l_mu_x += fabs(eq_plus + eq_min); // add surface charge
l_mu_x /= (2 * l_water); l_mu_x /= (2 * l_water);
eq_plus = eq_min = 0; eq_plus = eq_min = 0;
} }
@ -6140,6 +6155,8 @@ viscosity(cxxSurface *surf_ptr)
{ {
dw_t_visc = 0; dw_t_visc = 0;
t1 = l_moles / l_water; t1 = l_moles / l_water;
if (t1 < 1e-9)
continue;
l_z = fabs(s_x[i]->z); l_z = fabs(s_x[i]->z);
if (l_z) if (l_z)
f_z = (l_z * l_z + l_z) / 2; f_z = (l_z * l_z + l_z) / 2;
@ -6172,16 +6189,21 @@ viscosity(cxxSurface *surf_ptr)
if (!surf_ptr) s_x[i]->dw_t_visc = dw_t_visc + t3; if (!surf_ptr) s_x[i]->dw_t_visc = dw_t_visc + t3;
//output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc )); //output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc ));
} }
// parms for A... // parms for A and V_an. 7/26/24: added V_an calculation for gases z = 0
if ((l_z = s_x[i]->z) == 0) if ((l_z = s_x[i]->z) == 0)
continue;
Dw = s_x[i]->dw;
if (Dw)
{ {
Dw *= (0.89 / viscos_0 * tk_x / 298.15); if (s_x[i]->Jones_Dole[6])
if (s_x[i]->dw_t) {
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15); V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * l_moles;
m_an += l_moles;
}
continue;
} }
if ((Dw = s_x[i]->dw) == 0)
continue;
Dw *= (0.89 / viscos_0 * tk_x / 298.15);
if (s_x[i]->dw_t)
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15);
if (l_z < 0) if (l_z < 0)
{ {
if (!strcmp(s_x[i]->name, "Cl-")) if (!strcmp(s_x[i]->name, "Cl-"))
@ -6191,21 +6213,18 @@ viscosity(cxxSurface *surf_ptr)
V_an += V_Cl * l_moles; V_an += V_Cl * l_moles;
m_an += l_moles; m_an += l_moles;
} }
else// if (s_x[i]->Jones_Dole[6]) else if (s_x[i]->Jones_Dole[6])
{ {
V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * l_moles; V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * l_moles;
m_an += l_moles; m_an += l_moles;
} }
if (Dw) // anions for A...
{ m_min += l_moles;
// anions for A... t1 = l_moles * l_z;
m_min += l_moles; eq_min -= t1;
t1 = l_moles * l_z; eq_dw_min -= t1 / Dw;
eq_min -= t1;
eq_dw_min -= t1 / Dw;
}
} }
else if (Dw) else
{ {
// cations for A... // cations for A...
m_plus += l_moles; m_plus += l_moles;