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

First try at constant capacitance model.

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/trunk@9331 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2015-03-05 18:12:53 +00:00
parent 8914d2c58c
commit 6309425ae9
7 changed files with 142 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
basicsubs.cpp
View File

@ -540,7 +540,7 @@ diff_layer_total(const char *total_name, const char *surface_name)
int j;
for (j = 0; j < count_unknowns; j++)
{
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
if (x[j]->type != SURFACE_CB)
continue;
@ -582,7 +582,7 @@ diff_layer_total(const char *total_name, const char *surface_name)
*/
if (strcmp_nocase("psi", total_name) == 0)
{
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
return ((LDBLE) (x[j]->master[0]->s->la * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ));
@ -636,7 +636,7 @@ diff_layer_total(const char *total_name, const char *surface_name)
}
else if (strcmp_nocase("charge", total_name) == 0)
{
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL && dl_type_x == cxxSurface::NO_DL)
if ((use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM) && dl_type_x == cxxSurface::NO_DL)
{
return ((LDBLE) (x[j]->f));
}
@ -682,7 +682,7 @@ diff_layer_total(const char *total_name, const char *surface_name)
}
else if (strcmp_nocase("sigma", total_name) == 0)
{
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[j]->surface_charge);
if (dl_type_x != cxxSurface::NO_DL)

4
mainsubs.cpp
View File

@ -1679,7 +1679,7 @@ xsurface_save(int n_user)
}
comp_ptr->Set_charge_balance(charge);
}
else if (x[i]->type == SURFACE_CB && use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
else if (x[i]->type == SURFACE_CB && (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM))
{
cxxSurfaceCharge *charge_ptr = temp_surface.Find_charge(x[i]->surface_charge);
assert(charge_ptr);
@ -1977,7 +1977,7 @@ step_save_surf(int n_user)
/*
* Update grams
*/
if ((surface_ptr->Get_type() == cxxSurface::DDL || surface_ptr->Get_type() == cxxSurface::CD_MUSIC)
if ((surface_ptr->Get_type() == cxxSurface::DDL || surface_ptr->Get_type() == cxxSurface::CCM || surface_ptr->Get_type() == cxxSurface::CD_MUSIC)
&& surface_ptr->Get_related_rate() && use.Get_kinetics_ptr() != NULL)
{
for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++)

130
model.cpp
View File

@ -1981,35 +1981,61 @@ jacobian_sums(void)
*/
if (surface_unknown != NULL && dl_type_x == cxxSurface::NO_DL)
{
sinh_constant =
//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);
for (i = 0; i < count_unknowns; i++)
if (use.Get_surface_ptr()->Get_type() != cxxSurface::CCM)
{
cxxSurfaceCharge *charge_ptr = NULL;
if (x[i]->type == SURFACE_CB)
sinh_constant =
//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);
for (i = 0; i < count_unknowns; i++)
{
charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
}
if (x[i]->type == SURFACE_CB && charge_ptr->Get_grams() > 0)
{
for (j = 0; j < count_unknowns; j++)
cxxSurfaceCharge *charge_ptr = NULL;
if (x[i]->type == SURFACE_CB)
{
array[x[i]->number * (count_unknowns + 1) + j] *=
F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
}
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)
if (x[i]->type == SURFACE_CB && charge_ptr->Get_grams() > 0)
{
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);
for (j = 0; j < count_unknowns; j++)
{
array[x[i]->number * (count_unknowns + 1) + j] *=
F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
}
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) +
mu_unknown->number] -=
0.5 * sinh_constant / sqrt(mu_x) *
sinh(x[i]->master[0]->s->la * LOG_10);
}
}
}
}
else
{
for (i = 0; i < count_unknowns; i++)
{
cxxSurfaceCharge *charge_ptr = NULL;
if (x[i]->type == SURFACE_CB)
{
charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
}
if (x[i]->type == SURFACE_CB && charge_ptr->Get_grams() > 0)
{
for (j = 0; j < count_unknowns; j++)
{
array[x[i]->number * (count_unknowns + 1) + j] *=
F_C_MOL / (charge_ptr->Get_specific_area() *
charge_ptr->Get_grams());
}
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;
}
}
}
@ -3524,7 +3550,7 @@ reset(void)
/* recalculate g's for component */
if (dl_type_x != cxxSurface::NO_DL
&& (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL
&& (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM
|| (use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC
&& x[i]->type == SURFACE_CB2)))
{
@ -4349,6 +4375,60 @@ residuals(void)
iterations, x[i]->description, i, residual[i]));
converge = FALSE;
}
}
else if (x[i]->type == SURFACE_CB && use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
cxxSurfaceCharge *charge_ptr = use.Get_surface_ptr()->Find_charge(x[i]->surface_charge);
if (charge_ptr->Get_grams() == 0)
{
residual[i] = 0.0;
}
else if (dl_type_x != cxxSurface::NO_DL)
{
residual[i] = -x[i]->f;
}
else
{
residual[i] =
charge_ptr->Get_capacitance0() * x[i]->master[0]->s->la * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ -
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) */
)));
}
else
{
output_msg(sformatf( "\tResidual %e\n",
(double) x[i]->f));
}
output_msg(sformatf( "\t grams %g\n",
(double) charge_ptr->Get_grams()));
output_msg(sformatf( "\tCharge from potential %e eq\n",
(double) (charge_ptr->Get_capacitance0() * x[i]->master[0]->s->la * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ)));
output_msg(sformatf( "\t Psi %e\n",
(double) (x[i]->master[0]->s->la * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ)));
}
if (charge_ptr->Get_grams() > MIN_RELATED_SURFACE
&& fabs(residual[i]) > l_toler)
{
if (print_fail)
output_msg(sformatf(
"Failed Residual %d: %s %d %e\n", iterations,
x[i]->description, i, residual[i]));
converge = FALSE;
}
}
else if (x[i]->type == SURFACE_CB1)
{

6
prep.cpp
View File

@ -2836,7 +2836,7 @@ add_potential_factor(void)
error_msg(error_string, CONTINUE);
return(OK);
}
if (use.Get_surface_ptr()->Get_type() != cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() != cxxSurface::DDL && use.Get_surface_ptr()->Get_type() != cxxSurface::CCM)
return (OK);
sum_z = 0.0;
master_ptr = NULL;
@ -3069,7 +3069,7 @@ add_surface_charge_balance(void)
error_msg(error_string, CONTINUE);
return(OK);
}
if (use.Get_surface_ptr()->Get_type() != cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() != cxxSurface::DDL && use.Get_surface_ptr()->Get_type() != cxxSurface::CCM)
return (OK);
master_ptr = NULL;
/*
@ -3492,7 +3492,7 @@ setup_surface(void)
x[count_unknowns]->potential_unknown = NULL;
count_unknowns++;
/*if (use.Get_surface_ptr()->edl == FALSE) continue; */
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
/*
* Setup surface-potential unknown

17
print.cpp
View File

@ -1577,10 +1577,18 @@ print_surface(void)
*/
s_h2o->lm = s_h2o->la;
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
{
output_msg(sformatf("%-14s\n", "Diffuse Double Layer Surface-Complexation Model\n"));
}
else if (use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
output_msg(sformatf("%-14s\n", "Constant Capacitance Surface-Complexation Model\n"));
}
for (int j = 0; j < count_unknowns; j++)
{
/*if (use.Get_surface_ptr()->edl == TRUE) { */
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL)
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
if (x[j]->type != SURFACE_CB)
continue;
@ -1608,7 +1616,7 @@ print_surface(void)
(double) x[j]->f));
}
/*if (use.Get_surface_ptr()->edl == TRUE && diffuse_layer_x == FALSE) { */
if (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL && dl_type_x == cxxSurface::NO_DL)
if ((use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM) && dl_type_x == cxxSurface::NO_DL)
{
charge = x[j]->f;
}
@ -1641,6 +1649,11 @@ print_surface(void)
output_msg(sformatf("\tundefined sigma, C/m²\n"));
#endif
}
if (use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)
{
output_msg(sformatf("\t%11.3e capacitance, F/m^2\n",
(double) (charge_ptr->Get_capacitance0())));
}
output_msg(sformatf("\t%11.3e psi, V\n",
(double) (x[j]->master[0]->s->la * 2 * R_KJ_DEG_MOL *
tk_x * LOG_10 / F_KJ_V_EQ)));

15
read.cpp
View File

@ -7358,6 +7358,13 @@ read_surface(void)
case 13: /* ccm */
temp_surface.Set_type(cxxSurface::CCM);
copy_token(token1, &next_char);
if (charge_ptr == NULL)
{
error_msg("A surface must be defined before the capacitance is defined.\n",
CONTINUE);
input_error++;
break;
}
if (sscanf(token1.c_str(), SCANFORMAT, &dummy) != 1)
{
error_msg("Expected capacitance for constant_capacitance model.\n",
@ -7371,8 +7378,8 @@ read_surface(void)
}
/* constant capacitance model not implemented yet */
error_msg("Constant capacitance model not implemented.", CONTINUE);
input_error++;
//error_msg("Constant capacitance model not implemented.", CONTINUE);
//input_error++;
break;
case OPTION_DEFAULT:
@ -7623,7 +7630,7 @@ read_surface(void)
/*
* Make sure surface area is defined
*/
if (temp_surface.Get_type() == cxxSurface::DDL || temp_surface.Get_type() == cxxSurface::CD_MUSIC)
if (temp_surface.Get_type() == cxxSurface::DDL || temp_surface.Get_type() == cxxSurface::CCM || temp_surface.Get_type() == cxxSurface::CD_MUSIC)
{
for (size_t i = 0; i < temp_surface.Get_surface_charges().size(); i++)
{
@ -7657,7 +7664,7 @@ read_surface(void)
temp_surface.Set_dl_type(cxxSurface::NO_DL);
}
}
else if (temp_surface.Get_type() == cxxSurface::DDL)
else if (temp_surface.Get_type() == cxxSurface::DDL || temp_surface.Get_type() == cxxSurface::CCM)
{
/* all values of dl_type are valid */
}

4
step.cpp
View File

@ -553,12 +553,12 @@ add_surface(cxxSurface *surface_ptr)
}
}
}
if (surface_ptr->Get_type() != cxxSurface::DDL && surface_ptr->Get_type() != cxxSurface::CD_MUSIC)
if (surface_ptr->Get_type() != cxxSurface::DDL && surface_ptr->Get_type() != cxxSurface::CCM && surface_ptr->Get_type() != cxxSurface::CD_MUSIC)
return (OK);
for (size_t i = 0; i < surface_ptr->Get_surface_charges().size(); i++)
{
cxxSurfaceCharge *charge_ptr = &(surface_ptr->Get_surface_charges()[i]);
if (surface_ptr->Get_type() == cxxSurface::DDL || surface_ptr->Get_type() == cxxSurface::CD_MUSIC)
if (surface_ptr->Get_type() == cxxSurface::DDL || surface_ptr->Get_type() == cxxSurface::CCM || surface_ptr->Get_type() == cxxSurface::CD_MUSIC)
{
cb_x += charge_ptr->Get_charge_balance();
}