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

1088 lines
32 KiB
C++
Raw Blame History

#include "Phreeqc.h"
#include "phqalloc.h"
#include "Utils.h"
#include "Solution.h"
#define MAX_QUAD 20
#define K_POLY 5
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_all_g(void)
/* ---------------------------------------------------------------------- */
{
int converge, converge1;
LDBLE new_g, xd1;
LDBLE epsilon;
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* calculate g for each surface
*/
epsilon = convergence_tolerance;
if (convergence_tolerance >= 1e-8)
{
G_TOL = 1e-9;
}
else
{
G_TOL = 1e-10;
}
converge = TRUE;
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
continue;
if (debug_diffuse_layer == TRUE)
output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
std::map<LDBLE, cxxSurfDL> temp_g_map;
cxxSurfDL temp_g;
charge_ptr->Get_g_map()[0] = temp_g;
temp_g_map[0] = temp_g;
xd_global = exp(-2 * x[j]->master[0]->s->la * LOG_10);
/* alpha = 0.02935 @ 25; (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
/* 1000 J/kJ and 1000 L/m**3 */
//alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
// tk_x * 0.5);
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) * 1000.0 *
tk_x * 0.5);
/*
* calculate g for given surface for each species
*/
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
continue;
if (temp_g_map.find(s_x[i]->z) != temp_g_map.end())
continue;
z_global = s_x[i]->z;
if (charge_ptr->Get_grams() > 0.0)
{
if ((use.Get_surface_ptr()->Get_only_counter_ions() == false) ||
(((x[j]->master[0]->s->la > 0) && (z_global < 0))
|| ((x[j]->master[0]->s->la < 0) && (z_global > 0))))
{
if (xd_global > 0.1)
{
new_g = qromb_midpnt(charge_ptr, 1.0, xd_global);
}
else if (xd_global > 0.01)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, xd_global);
}
else if (xd_global > 0.001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, xd_global);
}
else if (xd_global > 0.0001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, xd_global);
}
else if (xd_global > 0.00001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
new_g += qromb_midpnt(charge_ptr, 0.0001, xd_global);
}
else if (xd_global > 0.000001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
new_g += qromb_midpnt(charge_ptr, 0.00001, xd_global);
}
else if (xd_global > 0.0000001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
new_g += qromb_midpnt(charge_ptr, 0.000001, xd_global);
}
else if (xd_global > 0.00000001)
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
new_g += qromb_midpnt(charge_ptr, 0.000001, .0000001);
new_g += qromb_midpnt(charge_ptr, 0.0000001, xd_global);
}
else
{
new_g = qromb_midpnt(charge_ptr, 1.0, 0.1);
new_g += qromb_midpnt(charge_ptr, 0.1, 0.01);
new_g += qromb_midpnt(charge_ptr, 0.01, .001);
new_g += qromb_midpnt(charge_ptr, 0.001, .0001);
new_g += qromb_midpnt(charge_ptr, 0.0001, .00001);
new_g += qromb_midpnt(charge_ptr, 0.00001, .000001);
new_g += qromb_midpnt(charge_ptr, 0.000001, .0000001);
new_g += qromb_midpnt(charge_ptr, 0.0000001, .00000001);
new_g += qromb_midpnt(charge_ptr, 0.00000001, xd_global);
}
}
else
{
new_g = 0;
}
}
else
{
new_g = 0.0;
}
if ((use.Get_surface_ptr()->Get_only_counter_ions()) && new_g < 0)
new_g = 0;
converge1 = TRUE;
if (fabs(new_g) >= 1.)
{
if (fabs((new_g - charge_ptr->Get_g_map()[z_global].Get_g()) / new_g) > epsilon)
{
converge1 = FALSE;
}
}
else
{
if (fabs(new_g - charge_ptr->Get_g_map()[z_global].Get_g()) > epsilon)
{
converge1 = FALSE;
}
}
if (converge1 == FALSE)
{
converge = FALSE;
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double) z_global,
(double) charge_ptr->Get_g_map()[z_global].Get_g(),
(double) new_g,
(double) (new_g - charge_ptr->Get_g_map()[z_global].Get_g())));
}
}
charge_ptr->Get_g_map()[z_global].Set_g(new_g);
if (new_g == 0)
{
charge_ptr->Get_g_map()[z_global].Set_dg(0.);
}
else
{
if (charge_ptr->Get_grams() > 0.0)
{
LDBLE dg = charge_ptr->Get_grams() *
charge_ptr->Get_specific_area() * alpha_global *
g_function(xd_global) / F_C_MOL;
dg *=
-2. / (exp(x[j]->master[0]->s->la * LOG_10) *
exp(x[j]->master[0]->s->la * LOG_10));
if ((xd_global - 1) < 0.0)
{
dg *= -1.0;
}
if (fabs(dg) < 1e-8)
{
xd1 = exp(-2 * 1e-3 * LOG_10);
new_g = qromb_midpnt(charge_ptr, 1.0, xd1);
dg = new_g / .001;
}
charge_ptr->Get_g_map()[z_global].Set_dg(dg);
}
else
{
charge_ptr->Get_g_map()[z_global].Set_dg(0.0);
}
}
temp_g_map[z_global] = charge_ptr->Get_g_map()[z_global];
}
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf("\nSurface component %d: charge,\tg,\tdg/dlny,\txd\n",
(int) charge_ptr->Get_g_map().size()));
std::map<LDBLE, cxxSurfDL>::iterator it;
for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
{
output_msg(sformatf(
"\t%12f\t%12.4e\t%12.4e\t%12.4e\n",
(double) it->first,
(double) it->second.Get_g(),
(double) it->second.Get_dg(),
(double) xd_global));
}
}
}
return (converge);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
g_function(LDBLE x_value)
/* ---------------------------------------------------------------------- */
{
LDBLE sum, return_value, sum1;
int i;
LDBLE ln_x_value;
if (equal(x_value, 1.0, G_TOL * 100) == TRUE)
return (0.0);
sum = 0.0;
ln_x_value = log(x_value);
cxxSurfaceCharge *charge_ptr = &(use.Get_surface_ptr()->Get_surface_charges()[0]);
std::map<LDBLE, cxxSurfDL>::iterator it = charge_ptr->Get_g_map().begin();
for ( ; it != charge_ptr->Get_g_map().end(); it++)
{
it->second.Set_psi_to_z(exp(ln_x_value * it->first) - 1.0);
}
for (i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
{
sum += s_x[i]->moles * charge_ptr->Get_g_map()[s_x[i]->z].Get_psi_to_z();
}
}
if (sum < 0.0)
{
sum = 0.0;
sum1 = 0.0;
output_msg(sformatf(
"Species\tmoles\tX**z-1\tsum\tsum charge\n"));
for (i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type < H2O && s_x[i]->z != 0.0)
{
sum += s_x[i]->moles * (pow(x_value, s_x[i]->z) - 1.0);
sum1 += s_x[i]->moles * s_x[i]->z;
output_msg(sformatf( "%s\t%e\t%e\t%e\t%e\n",
s_x[i]->name, (double) s_x[i]->moles,
(double) (pow((LDBLE) x_value, (LDBLE) s_x[i]->z) -
1.0), (double) sum, (double) sum1));
}
}
error_string = sformatf( "Negative sum in g_function, %e\t%e.",
(double) sum, (double) x_value);
error_msg(error_string, CONTINUE);
error_string = sformatf(
"Solutions must be charge balanced, charge imbalance is %e\n",
(double) sum1);
error_msg(error_string, STOP);
}
return_value =
(exp(ln_x_value * z_global) -
1) / sqrt((x_value * x_value * mass_water_aq_x * sum));
return (return_value);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
polint(LDBLE * xa, LDBLE * ya, int n, LDBLE xv, LDBLE * yv, LDBLE * dy)
/* ---------------------------------------------------------------------- */
{
int i, m, ns;
LDBLE den, dif, dift, ho, hp, w;
LDBLE *c, *d;
ns = 1;
dif = fabs(xv - xa[1]);
/*
* Malloc work space
*/
c = (LDBLE *) PHRQ_malloc((size_t) (n + 1) * sizeof(LDBLE));
if (c == NULL)
malloc_error();
d = (LDBLE *) PHRQ_malloc((size_t) (n + 1) * sizeof(LDBLE));
if (d == NULL)
malloc_error();
for (i = 1; i <= n; i++)
{
dift = fabs(xv - xa[i]);
if (dift < dif)
{
ns = i;
dif = dift;
}
c[i] = ya[i];
d[i] = ya[i];
}
*yv = ya[ns--];
for (m = 1; m < n; m++)
{
for (i = 1; i <= n - m; i++)
{
ho = xa[i] - xv;
hp = xa[i + m] - xv;
w = c[i + 1] - d[i];
if ((den = ho - hp) == 0.0)
{
error_msg("In subroutine polint.", STOP);
}
den = w / den;
d[i] = hp * den;
c[i] = ho * den;
}
if (2 * ns < (n - m))
{
*dy = c[ns + 1];
}
else
{
*dy = d[ns--];
}
*yv += *dy;
/* *yv += (*dy = (2 * ns < (n-m) ? c[ns+1] : d[ns--])); */
}
c = (LDBLE *) free_check_null(c);
d = (LDBLE *) free_check_null(d);
return;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
midpnt(LDBLE x1, LDBLE x2, int n)
/* ---------------------------------------------------------------------- */
{
LDBLE xv, tnm, sum, del, ddel;
int it, j;
if (n == 1)
{
midpoint_sv = (x2 - x1) * g_function(0.5 * (x1 + x2));
return (midpoint_sv);
}
else
{
for (it = 1, j = 1; j < n - 1; j++)
it *= 3;
tnm = (LDBLE) it;
del = (x2 - x1) / (3 * tnm);
ddel = del + del;
xv = x1 + 0.5 * del;
sum = 0.0;
for (j = 1; j <= it; j++)
{
#if defined(PHREEQCI_GUI)
PhreeqcIWait(this);
#endif
sum += g_function(xv);
xv += ddel;
sum += g_function(xv);
xv += del;
}
midpoint_sv = (midpoint_sv + (x2 - x1) * sum / tnm) / 3.0;
return midpoint_sv;
}
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
qromb_midpnt(cxxSurfaceCharge *charge_ptr, LDBLE x1, LDBLE x2)
/* ---------------------------------------------------------------------- */
{
LDBLE ss, dss;
LDBLE sv[MAX_QUAD + 2], h[MAX_QUAD + 2];
int j;
h[0] = 1.0;
sv[0] = midpnt(x1, x2, 1);
for (j = 1; j < MAX_QUAD; j++)
{
sv[j] = midpnt(x1, x2, j + 1);
h[j] = h[j - 1] / 9.0;
if (fabs(sv[j] - sv[j - 1]) <= G_TOL * fabs(sv[j]))
{
sv[j] *= charge_ptr->Get_grams() * charge_ptr->Get_specific_area() * alpha_global / F_C_MOL; /* (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
if ((x2 - 1) < 0.0)
sv[j] *= -1.0;
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"Iterations in qromb_midpnt: %d\n", j));
}
return (sv[j]);
}
if (j >= K_POLY - 1)
{
polint(&h[j - K_POLY], &sv[j - K_POLY], K_POLY, 0.0, &ss, &dss);
if (fabs(dss) <= G_TOL * fabs(ss) || fabs(dss) < G_TOL)
{
ss *= charge_ptr->Get_grams() * charge_ptr->Get_specific_area() * alpha_global / F_C_MOL; /* (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
if ((x2 - 1) < 0.0)
ss *= -1.0;
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"Iterations in qromb_midpnt: %d\n", j));
}
return (ss);
}
}
}
error_string = sformatf(
"\nToo many iterations integrating diffuse layer.\n");
error_msg(error_string, STOP);
return (-999.9);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_init_g(void)
/* ---------------------------------------------------------------------- */
{
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* calculate g for each surface
*/
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
continue;
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
xd_global = exp(-2 * x[j]->master[0]->s->la * LOG_10);
/* alpha = 0.02935 @ 25; (ee0RT/2)**1/2, (L/mol)**1/2 C / m**2 */
/* second 1000 is liters/m**3 */
//alpha_global = sqrt(EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// 1000.0 * tk_x * 0.5);
alpha_global = sqrt(eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
1000.0 * tk_x * 0.5);
if (charge_ptr->Get_g_map().size() == 0)
{
cxxSurfDL temp_g;
charge_ptr->Get_g_map()[0.0] = temp_g;
}
/*
* calculate g for given surface for each species
*/
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
continue;
if (charge_ptr->Get_g_map().find(s_x[i]->z) == charge_ptr->Get_g_map().end())
{
cxxSurfDL temp_g;
/* save g for charge */
if (charge_ptr->Get_grams() > 0.0)
{
temp_g.Set_g(2 * alpha_global * sqrt(mu_x) * (pow(xd_global, s_x[i]->z / 2.0) - 1) *
charge_ptr->Get_grams() *
charge_ptr->Get_specific_area() / F_C_MOL);
temp_g.Set_dg(-s_x[i]->z);
if (use.Get_surface_ptr()->Get_only_counter_ions() &&
temp_g.Get_g() < 0)
{
temp_g.Set_g(0);
temp_g.Set_dg(0);
}
}
else
{
temp_g.Set_g(0);
temp_g.Set_dg(-s_x[i]->z);
}
charge_ptr->Get_g_map()[s_x[i]->z] = temp_g;
}
{
int is = s_x[i]->number;
assert (is < (int) s_diff_layer.size());
// species found in diff_layer
s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0);
s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0);
}
}
if (debug_diffuse_layer == TRUE)
{
output_msg(sformatf(
"\nSurface component %d: charge,\tg,\tdg\n",
(int) charge_ptr->Get_g_map().size()));
std::map<LDBLE, cxxSurfDL>::iterator it;
for (it = charge_ptr->Get_g_map().begin(); it != charge_ptr->Get_g_map().end(); it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) it->first,
(double) it->second.Get_g(),
(double) it->second.Get_dg()));
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
initial_surface_water(void)
/* ---------------------------------------------------------------------- */
{
/*
* In initial surface calculation, need to calculate
* mass of water in diffuse layer.
* diffuse layer water + aqueous solution water = bulk water.
* Ionic strength is fixed, so diffuse-layer water will not change
*/
LDBLE debye_length, b, r, rd, ddl_limit, rd_limit, fraction, sum_surfs, l_s;
LDBLE damp_aq;
/*
* Debye length = 1/k = sqrt[eta*eta_zero*R*T/(2*F**2*mu_x*1000)], Dzombak and Morel, p 36
*
* 1000 converts kJ to J; 1000 converts Liters to meter**3; debye_length is in meters.
*/
//debye_length = (EPSILON * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
// / (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
debye_length = (eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * tk_x)
/ (2. * F_C_MOL * F_C_MOL * mu_x * 1000.);
debye_length = sqrt(debye_length);
/* ddl is at most the fraction ddl_limit of bulk water */
ddl_limit = use.Get_surface_ptr()->Get_DDL_limit();
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
if (use.Get_surface_ptr()->Get_debye_lengths() > 0)
{
sum_surfs = 0.0;
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);
sum_surfs +=
charge_ptr->Get_specific_area() *
charge_ptr->Get_grams();
}
rd = debye_length * use.Get_surface_ptr()->Get_debye_lengths();
use.Get_surface_ptr()->Set_thickness(rd);
if (!sum_surfs)
{
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);
charge_ptr->Set_mass_water(0);
}
}
else if (state == INITIAL_SURFACE)
{
/* distribute water over DDL (rd) and free pore (r - rd) */
/* find r: free pore (m3) = pi * (r - rd)^2 * L, where L = A / (2*pi*r),
A = sum_surfs = sum of the surface areas */
b = -2 * (rd + use.Get_solution_ptr()->Get_mass_water() / (1000 * sum_surfs));
r = 0.5 * (-b + sqrt(b * b - 4 * rd * rd));
/* DDL (m3) = pi * (r^2 - (r - rd)^2) * L */
rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
/* rd should be smaller than r and the ddl limit */
if (rd > rd_limit)
{
mass_water_surfaces_x =
use.Get_solution_ptr()->Get_mass_water() * ddl_limit / (1 - ddl_limit);
r = 0.002 * (mass_water_surfaces_x +
use.Get_solution_ptr()->Get_mass_water()) / sum_surfs;
rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
rd = rd_limit;
use.Get_surface_ptr()->Set_thickness(rd);
}
else
mass_water_surfaces_x =
(r * r / pow(r - rd, 2) - 1) * use.Get_solution_ptr()->Get_mass_water();
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();
charge_ptr->Set_mass_water(mass_water_surfaces_x * l_s / sum_surfs);
}
}
else
{
r = 0.002 * mass_water_bulk_x / sum_surfs;
rd_limit = (1 - sqrt(1 - ddl_limit)) * r;
if (rd > rd_limit)
{
rd = rd_limit;
use.Get_surface_ptr()->Set_thickness(rd);
fraction = ddl_limit;
}
else
fraction = 1 - pow(r - rd, 2) / (r * r);
damp_aq = 1.0;
if (g_iterations > 10)
damp_aq = 0.2;
else if (g_iterations > 5)
damp_aq = 0.5;
mass_water_surfaces_x = damp_aq * fraction * mass_water_bulk_x +
(1 - damp_aq) * mass_water_surfaces_x;
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();
charge_ptr->Set_mass_water(mass_water_surfaces_x * l_s / sum_surfs);
}
}
}
else
{
/* take constant thickness of, default 1e-8 m (100 Angstroms) */
mass_water_surfaces_x = 0.0;
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);
charge_ptr->Set_mass_water(charge_ptr->Get_specific_area() *
charge_ptr->Get_grams() * use.Get_surface_ptr()->Get_thickness() *
1000);
mass_water_surfaces_x += charge_ptr->Get_mass_water();
}
}
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC)
mass_water_bulk_x = mass_water_aq_x + mass_water_surfaces_x;
else
{
/* for variable distribution of water over DDL and bulk... */
if (state > INITIAL_SURFACE)
mass_water_aq_x = mass_water_bulk_x - mass_water_surfaces_x;
else
mass_water_bulk_x = mass_water_aq_x + mass_water_surfaces_x;
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
sum_diffuse_layer(cxxSurfaceCharge *charge_ptr)
/* ---------------------------------------------------------------------- */
{
LDBLE mass_water_surface;
LDBLE molality, moles_excess, moles_surface;
if (use.Get_surface_ptr() == NULL)
return (OK);
/*
* Find position of component in list of components
*/
/*
* Loop through all surface components, calculate each H2O surface (diffuse layer),
* H2O aq, and H2O bulk (diffuse layers plus aqueous).
*/
count_elts = 0;
paren_count = 0;
mass_water_surface = charge_ptr->Get_mass_water();
for (int j = 0; j < (int)this->s_x.size(); j++)
{
if (s_x[j]->type > HPLUS)
continue;
molality = under(s_x[j]->lm);
LDBLE g = charge_ptr->Get_g_map()[s_x[j]->z].Get_g();
if (s_x[j]->erm_ddl != 1)
{
LDBLE ratio_aq = mass_water_surface / mass_water_aq_x;
LDBLE g2 = g / ratio_aq + 1;
g = ratio_aq * (g2 * s_x[j]->erm_ddl - 1);
}
moles_excess = mass_water_aq_x * molality * g;
moles_surface = mass_water_surface * molality + moles_excess;
/*
* Accumulate elements in diffuse layer
*/
add_elt_list(s_x[j]->next_elt, moles_surface);
}
add_elt_list(s_h2o->next_elt, mass_water_surface / gfw_water);
elt_list_combine();
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_all_donnan(void)
/* ---------------------------------------------------------------------- */
{
bool converge;
int cd_m;
LDBLE new_g, f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq, ratio_aq_tot;
LDBLE new_g2, f_psi2, surf_chrg_eq2, psi_avg2, dif, var1;
if (use.Get_surface_ptr() == NULL)
return (OK);
//f_sinh = sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// tk_x * mu_x);
f_sinh = sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
tk_x * mu_x);
/*
* calculate g for each surface...
*/
if (!calculating_deriv || use.Get_surface_ptr()->Get_debye_lengths()) // DL_pitz
initial_surface_water();
converge = TRUE;
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
continue;
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf( "Calc_all_g, X[%d]\n", j));
/*
* sum eq of each charge number in solution...
*/
std::map<LDBLE, LDBLE>::iterator it;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
it->second = 0.0;
}
charge_group_map.clear();
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
continue;
charge_group_map[s_x[i]->z] += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
}
/* find surface charge from potential... */
A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC)
{
f_psi = x[j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2;
cd_m = 1;
} else
{
f_psi = x[j]->master[0]->s->la * LOG_10;
cd_m = -1;
}
surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
if (fabs(surf_chrg_eq) > 5e3)
{
surf_chrg_eq = (surf_chrg_eq < 0 ? -5e3 : 5e3);
var1 = surf_chrg_eq / (A_surf * f_sinh / F_C_MOL);
var1 = (var1 + sqrt(var1 * var1 + 1));
f_psi = (var1 > 1e-8 ? log(var1) : -18.4);
surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
x[j]->master[0]->s->la = f_psi / LOG_10;
}
/* also for the derivative... */
dif = 1e-5;
f_psi2 = f_psi + dif;
surf_chrg_eq2 = A_surf * f_sinh * sinh(f_psi2) / F_C_MOL;
/* find psi_avg that matches surface charge... */
psi_avg = calc_psi_avg(charge_ptr, surf_chrg_eq);
psi_avg2 = calc_psi_avg(charge_ptr, surf_chrg_eq2);
/*output_msg(sformatf( "psi's %e %e %e\n", f_psi, psi_avg, surf_chrg_eq)); */
/* fill in g's */
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
ratio_aq_tot = charge_ptr->Get_mass_water() / mass_water_bulk_x;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
LDBLE z = it->first;
if (!ratio_aq)
{
charge_ptr->Get_g_map()[z].Set_g(0);
charge_ptr->Get_g_map()[z].Set_dg(0);
charge_ptr->Get_z_gMCD_map()[z] = 0;
converge = true;
continue;
}
new_g = ratio_aq * (exp(cd_m * z * psi_avg) - 1);
if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
//((surf_chrg_eq < 0 && z < 0)
// || (surf_chrg_eq > 0 && z > 0)))
new_g = -ratio_aq;
if (new_g <= -ratio_aq)
new_g = -ratio_aq + G_TOL * 1e-3;
new_g2 = ratio_aq * (exp(cd_m * z * psi_avg2) - 1);
if (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0)
//((surf_chrg_eq < 0 && z < 0)
// || (surf_chrg_eq > 0 && z > 0)))
new_g2 = -ratio_aq;
if (new_g2 <= -ratio_aq)
new_g2 = -ratio_aq + G_TOL * 1e-3;
if (fabs(new_g) >= 1)
{
if (fabs((new_g - charge_ptr->Get_g_map()[z].Get_g()) / new_g) > convergence_tolerance)
{
converge = FALSE;
}
}
else
{
if (fabs(new_g - charge_ptr->Get_g_map()[z].Get_g()) > convergence_tolerance)
{
converge = FALSE;
}
}
charge_ptr->Get_g_map()[z].Set_g(new_g);
if (new_g != 0)
{
charge_ptr->Get_g_map()[z].Set_dg((new_g2 - new_g) / dif);
}
else
{
charge_ptr->Get_g_map()[z].Set_dg(-z);
}
/* save Boltzmann factor * water fraction for MCD calc's in transport */
if (converge)
{
if (use.Get_surface_ptr()->Get_only_counter_ions())
{
if (surf_chrg_eq * z > 0) // co-ions are not in the DL
charge_ptr->Get_z_gMCD_map()[z] = 0;
else // assume that counter-ions have the free water conc for diffusion
charge_ptr->Get_z_gMCD_map()[z] = ratio_aq_tot;
}
else
charge_ptr->Get_z_gMCD_map()[z] = (new_g / ratio_aq + 1) * ratio_aq_tot;
}
}
if (debug_diffuse_layer == TRUE)
{
std::string name = x[j]->master[0]->elt->name;
Utilities::replace("_psi", "", name);
output_msg(sformatf(
"\nDonnan all on %s (%d): charge, \tg, \tdg, Psi_surface = %8f V. \n",
name.c_str(), (int) charge_ptr->Get_g_map().size(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
tk_x / F_KJ_V_EQ));
for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
i_it != charge_ptr->Get_g_map().end(); i_it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) i_it->first,
(double) i_it->second.Get_g(),
(double) i_it->second.Get_dg()));
}
}
}
return (converge);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
calc_init_donnan(void)
/* ---------------------------------------------------------------------- */
{
LDBLE f_psi, surf_chrg_eq, psi_avg, f_sinh, A_surf, ratio_aq;
if (use.Get_surface_ptr() == NULL)
return (OK);
f_sinh =
//sqrt(8000.0 * EPSILON * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
// tk_x * mu_x);
sqrt(8000.0 * eps_r * EPSILON_ZERO * (R_KJ_DEG_MOL * 1000.0) *
tk_x * mu_x);
if (convergence_tolerance >= 1e-8)
{
G_TOL = 1e-9;
}
else
{
G_TOL = 1e-13;
}
/*
* sum eq of each charge number in solution...
*/
charge_group_map.clear();
charge_group_map[0.0] = 0.0;
for (int i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type > HPLUS)
continue;
if (charge_group_map.find(s_x[i]->z) != charge_group_map.end())
{
charge_group_map.find(s_x[i]->z)->second += s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
}
else
{
charge_group_map[s_x[i]->z] = s_x[i]->z * s_x[i]->moles * s_x[i]->erm_ddl;
}
}
/*
* calculate g for each surface...
*/
for (int j = 0; j < count_unknowns; j++)
{
if (x[j]->type != SURFACE_CB)
continue;
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
charge_ptr->Get_g_map().clear();
/* find surface charge from potential... */
A_surf = charge_ptr->Get_specific_area() * charge_ptr->Get_grams();
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC)
{
f_psi = x[j + 2]->master[0]->s->la * LOG_10; /* -FPsi/RT */
f_psi = f_psi / 2;
} else
f_psi = x[j]->master[0]->s->la * LOG_10;
surf_chrg_eq = A_surf * f_sinh * sinh(f_psi) / F_C_MOL;
/* find psi_avg that matches surface charge... */
/* psi_avg = calc_psi_avg (0);
appt 7/9/8... may have to change above one */
psi_avg = calc_psi_avg(charge_ptr, 0 * surf_chrg_eq);
/* fill in g's */
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
std::map<LDBLE, LDBLE>::iterator kit;
for (kit = charge_group_map.begin(); kit != charge_group_map.end(); kit++)
{
LDBLE z = kit->first;
LDBLE eq = kit->second;
charge_ptr->Get_g_map()[z].Set_g(ratio_aq * (exp(-z * psi_avg) - 1));
if (use.Get_surface_ptr()->Get_only_counter_ions()
&& ((surf_chrg_eq < 0 && z < 0)
|| (surf_chrg_eq > 0 && z > 0)))
charge_ptr->Get_g_map()[z].Set_g(-ratio_aq);
if (charge_ptr->Get_g_map()[z].Get_g() != 0)
{
charge_ptr->Get_g_map()[z].Set_dg(-A_surf * f_sinh * cosh(f_psi) /
(eq * F_C_MOL));
}
else
{
charge_ptr->Get_g_map()[z].Set_dg(-z);
}
/* save g for species */
for (int i = 0; i < (int)this->s_x.size(); i++)
{
int is = s_x[i]->number;
assert (is < (int) s_diff_layer.size());
s_diff_layer[is][charge_ptr->Get_name()].Set_g_moles(0.0);
s_diff_layer[is][charge_ptr->Get_name()].Set_dg_g_moles(0.0);
}
}
if (debug_diffuse_layer == TRUE)
{
std::string name = x[j]->master[0]->elt->name;
Utilities::replace("_psi", "", name);
output_msg(sformatf(
"\nDonnan init on %s : charge, \tg, \tdg, Psi_surface = %8f V. \n",
name.c_str(),
x[j]->master[0]->s->la * 2 * LOG_10 * R_KJ_DEG_MOL *
tk_x / F_KJ_V_EQ));
for (std::map<LDBLE, cxxSurfDL>::iterator i_it = charge_ptr->Get_g_map().begin();
i_it != charge_ptr->Get_g_map().end(); i_it++)
{
output_msg(sformatf( "\t%12f\t%12.4e\t%12.4e\n",
(double) i_it->first,
(double) i_it->second.Get_g(),
(double) i_it->second.Get_dg()));
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
/* ---------------------------------------------------------------------- */
{
/*
* calculate the average (F * Psi / RT) that lets the DL charge counter the surface charge
*/
LDBLE fd, fd1, p, temp, ratio_aq;
ratio_aq = charge_ptr->Get_mass_water() / mass_water_aq_x;
p = 0;
if (surf_chrg_eq == 0 || ratio_aq == 0)
return (0.0);
else if (surf_chrg_eq < 0)
p = -0.5 * log(-surf_chrg_eq * ratio_aq / mu_x + 1);
else if (surf_chrg_eq > 0)
p = 0.5 * log(surf_chrg_eq * ratio_aq / mu_x + 1);
/*
* Optimize p in SS{s_x[i]->moles * z_i * g(p)} = -surf_chrg_eq
* g(p) = exp(-p * z_i) * ratio_aq
* Elsewhere in PHREEQC, g is the excess, after subtraction of conc's for p = 0:
* g(p) = (exp(-p *z_i) - 1) * ratio_aq
*/
int l_iter = 0;
do
{
fd = surf_chrg_eq;
fd1 = 0.0;
std::map<LDBLE, LDBLE>::iterator it;
for (it = charge_group_map.begin(); it != charge_group_map.end(); it++)
{
LDBLE z = it->first;
if (!z || (use.Get_surface_ptr()->Get_only_counter_ions() && surf_chrg_eq * z > 0))
continue;
LDBLE eq = it->second;
/* multiply with ratio_aq for multiplier options cp and cm
in calc_all_donnan (not used now)... */
temp = exp(-z * p) * ratio_aq;
fd += eq * temp;
fd1 -= z * eq * temp;
}
fd /= -fd1;
p += (fd > 1) ? 1 : ((fd < -1) ? -1 : fd);
if (fabs(p) < G_TOL)
p = 0.0;
l_iter++;
if (l_iter > 50)
{
pr.all = TRUE;
pr.exchange = TRUE;
pr.headings = TRUE;
pr.pp_assemblage = TRUE;
pr.species = TRUE;
pr.surface = TRUE;
pr.totals = TRUE;
print_all();
error_string = sformatf(
"\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n",
(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
error_msg(error_string, STOP);
}
}
while (fabs(fd) > 1e-12 && p != 0.0);
if (debug_diffuse_layer == TRUE)
output_msg(sformatf(
"iter in calc_psi_avg = %d. g(+1) = %8f. surface charge = %12.4e.\n",
l_iter, (double) (exp(-p) - 1), (double) surf_chrg_eq));
return (p);
}
Reference in New Issue View Git Blame Copy Permalink