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

Merge commit 'd7188d25f30f33d2158082341224dda94f4bc1cc'

Browse Source
This commit is contained in:
Scott R Charlton 2019-11-19 09:40:14 -07:00
commit 9a6afb37c4
3 changed files with 98 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
phreeqcpp/cvode.cpp
View File

@ -3590,7 +3590,8 @@ CVsldet(CVodeMem cv_mem)
realtype smink, smaxk, sumrat, sumrsq, vmin, vmax, drrmax, adrr;
realtype /*small_cvode,*/ tem, sqmax, saqk, qp, s, sqmaxk, saqj, sqmin;
realtype rsa, rsb, rsc, rsd, rse, rd1a, rd1b, rd1c, rd1d;
realtype rd2a, rd2b, rd2c, rd3a, rd3b, cest1, corr1;
//realtype rd2a, rd2b, rd2c, rd3a, rd3b, cest1, corr1;
realtype rd2a, rd2b, rd2c, rd3a, cest1, corr1;
realtype ratp, ratm, qfac1, qfac2, bb, rrb;
/* The following are cutoffs and tolerances used by this routine */
@ -3854,7 +3855,7 @@ CVsldet(CVodeMem cv_mem)
rd2b = rd1b - rd1c;
rd2c = rd1c - rd1d;
rd3a = rd2a - rd2b;
rd3b = rd2b - rd2c;
//rd3b = rd2b - rd2c;
/* rd3b = rd3b; */
if (ABS(rd1b) < TINY * smax[k])

163
phreeqcpp/model.cpp
View File

@ -687,98 +687,99 @@ gammas(LDBLE mu)
*/
if (calculating_deriv)
continue;
LDBLE coef, z;
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
if (s_x[i]->rxn_x->token[j].s->type == EX)
LDBLE coef = 0, z = 0;
for (j = 1; s_x[i]->rxn_x->token[j].s != NULL; j++)
{
s_x[i]->alk =
s_x[i]->rxn_x->token[j].s->primary->unknown->moles;
//break;
if (s_x[i]->rxn_x->token[j].s->type == EX)
{
s_x[i]->alk =
s_x[i]->rxn_x->token[j].s->primary->unknown->moles;
//break;
}
else if (s_x[i]->rxn_x->token[j].s->type <= HPLUS)
{
coef = s_x[i]->rxn_x->token[j].coef;
z = s_x[i]->rxn_x->token[j].s->z;
}
}
else if (s_x[i]->rxn_x->token[j].s->type <= HPLUS)
if (!use.Get_exchange_ptr()->Get_pitzer_exchange_gammas())
{
coef = s_x[i]->rxn_x->token[j].coef;
z = s_x[i]->rxn_x->token[j].s->z;
}
}
if (!use.Get_exchange_ptr()->Get_pitzer_exchange_gammas())
{
if (s_x[i]->primary != NULL)
{
s_x[i]->lg = 0.0;
s_x[i]->dg = 0.0;
}
else
{
if (s_x[i]->alk <= 0)
if (s_x[i]->primary != NULL)
{
s_x[i]->lg = 0.0;
s_x[i]->dg = 0.0;
}
else
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = 0.0;
{
if (s_x[i]->alk <= 0)
s_x[i]->lg = 0.0;
else
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = 0.0;
}
}
}
else if (s_x[i]->exch_gflag == 1 && s_x[i]->alk > 0)
{
/* Davies */
s_x[i]->lg = -coef * z * z * a *
(muhalf / (1.0 + muhalf) - 0.3 * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg =
c1 * coef * z * z * s_x[i]->moles;
}
else if (s_x[i]->exch_gflag == 2 && s_x[i]->alk > 0)
{
/* Extended D-H, WATEQ D-H */
s_x[i]->lg = coef * (-a * muhalf * z * z /
(1.0 + s_x[i]->dha * b * muhalf) + s_x[i]->dhb * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = coef * (c2 * z * z /
((1.0 + s_x[i]->dha * b * muhalf) * (1.0 + s_x[i]->dha * b * muhalf)) +
s_x[i]->dhb) * LOG_10 * s_x[i]->moles;
}
else if (s_x[i]->exch_gflag == 7 && s_x[i]->alk > 0)
{
if (llnl_count_temp > 0)
else if (s_x[i]->exch_gflag == 1 && s_x[i]->alk > 0)
{
s_x[i]->lg =
coef * (-a_llnl * muhalf * z * z /
(1.0 + s_x[i]->dha * b_llnl * muhalf) +
bdot_llnl * mu) +
/* Davies */
s_x[i]->lg = -coef * z * z * a *
(muhalf / (1.0 + muhalf) - 0.3 * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg =
coef * (c2_llnl * z * z /
((1.0 + s_x[i]->dha * b_llnl * muhalf) * (1.0 + s_x[i]->dha * b_llnl * muhalf)) +
bdot_llnl) * LOG_10 * s_x[i]->moles;
c1 * coef * z * z * s_x[i]->moles;
}
else
else if (s_x[i]->exch_gflag == 2 && s_x[i]->alk > 0)
{
error_msg("LLNL_AQUEOUS_MODEL_PARAMETERS not defined.",
STOP);
/* Extended D-H, WATEQ D-H */
s_x[i]->lg = coef * (-a * muhalf * z * z /
(1.0 + s_x[i]->dha * b * muhalf) + s_x[i]->dhb * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = coef * (c2 * z * z /
((1.0 + s_x[i]->dha * b * muhalf) * (1.0 + s_x[i]->dha * b * muhalf)) +
s_x[i]->dhb) * LOG_10 * s_x[i]->moles;
}
}
else
{
/*
* Master species is a dummy variable with meaningless activity and mass
*/
if (s_x[i]->primary != NULL)
else if (s_x[i]->exch_gflag == 7 && s_x[i]->alk > 0)
{
s_x[i]->lg = 0.0;
s_x[i]->dg = 0.0;
}
else
{
if (s_x[i]->alk <= 0)
s_x[i]->lg = 0.0;
if (llnl_count_temp > 0)
{
s_x[i]->lg =
coef * (-a_llnl * muhalf * z * z /
(1.0 + s_x[i]->dha * b_llnl * muhalf) +
bdot_llnl * mu) +
log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg =
coef * (c2_llnl * z * z /
((1.0 + s_x[i]->dha * b_llnl * muhalf) * (1.0 + s_x[i]->dha * b_llnl * muhalf)) +
bdot_llnl) * LOG_10 * s_x[i]->moles;
}
else
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = 0.0;
{
error_msg("LLNL_AQUEOUS_MODEL_PARAMETERS not defined.",
STOP);
}
}
else
{
/*
* Master species is a dummy variable with meaningless activity and mass
*/
if (s_x[i]->primary != NULL)
{
s_x[i]->lg = 0.0;
s_x[i]->dg = 0.0;
}
else
{
if (s_x[i]->alk <= 0)
s_x[i]->lg = 0.0;
else
s_x[i]->lg = log10(fabs(s_x[i]->equiv) / s_x[i]->alk);
s_x[i]->dg = 0.0;
}
}
if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)
gammas_a_f(i); // appt
}
if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)
gammas_a_f(i); // appt
break;
case 5: /* Always 1.0 */
s_x[i]->lg = 0.0;
@ -879,7 +880,8 @@ int Phreeqc::gammas_a_f(int i1)
/* ------------------------------------------------------------------------------- */
{
int i, j;
LDBLE d2, d3, coef = 0, sum = 0;
//LDBLE d2, d3, coef = 0, sum = 0;
LDBLE d2, d3, sum = 0;
char name[MAX_LENGTH];
//struct master *m_ptr;
@ -2219,7 +2221,7 @@ mb_ss(void)
{
cxxSS *ss_ptr = ss_ptrs[i];
total_moles = 0;
bool ss_in = true;
//bool ss_in = true;
for (size_t j = 0; j < ss_ptr->Get_ss_comps().size(); j++)
{
int l;
@ -3020,7 +3022,8 @@ ss_binary(cxxSS *ss_ptr)
/* ---------------------------------------------------------------------- */
{
LDBLE nb, nc, n_tot, xb, xc, dnb, dnc, l_a0, l_a1;
LDBLE xb2, xb3, xb4, xc2, xc3;
//LDBLE xb2, xb3, xb4, xc2, xc3;
LDBLE xb2, xb3, xc2;
LDBLE xb1, xc1;
/*
* component 0 is major component
@ -3107,10 +3110,10 @@ ss_binary(cxxSS *ss_ptr)
comp1_ptr->Get_log10_lambda();
xc2 = xc * xc;
xc3 = xc2 * xc;
//xc3 = xc2 * xc;
xb2 = xb * xb;
xb3 = xb2 * xb;
xb4 = xb3 * xb;
//xb4 = xb3 * xb;
/* xb4 = xb4; */
/* xc3 = xc3; */

19
phreeqcpp/transport.cpp
View File

@ -2025,7 +2025,7 @@ fill_spec(int l_cell_no, int ref_cell)
malloc_error();
sol_D[l_cell_no].spec_size = i3 + count_spec + 1 + size_xt;
}
for (i1; i1 < i2; i1++)
for (; i1 < i2; i1++) // i1 is loop variable
{
memmove(&sol_D[l_cell_no].spec[i1], &sol_D[ref_cell].spec[i1], sizeof(struct spec));
sol_D[l_cell_no].spec[i1].c = 0.0;
@ -2098,7 +2098,7 @@ fill_spec(int l_cell_no, int ref_cell)
sol_D[i1].spec_size = i2 + size_xt;
}
i2--;
for (i2; i2 > i3; i2--)
for (; i2 > i3; i2--) // i2 is loop variable
sol_D[i1].spec[i2] = sol_D[i1].spec[i2 - 1];
memmove(&sol_D[i1].spec[i2], &sol_D[l_cell_no].spec[i2], sizeof(struct spec));
@ -2189,9 +2189,9 @@ diffuse_implicit(LDBLE DDt, int stagnant)
int i, icell, cp, comp;
// ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells);
int ifirst, ilast;
int i_1, i0, i1, i2;
int i_1, i0, i1, i2 = 0;
double mfr, mfr1, max_b = 0, b, grad, dVc, j_0e, min_dif_M = pow(10, min_dif_LM);
LDBLE dum1, dum2, dum_stag, min_mol;
LDBLE dum1, dum2, dum_stag = 0.0, min_mol;
LDBLE dum = 0;
cxxSurfaceCharge * charge_ptr = NULL;
@ -2789,10 +2789,15 @@ diffuse_implicit(LDBLE DDt, int stagnant)
for (i = ifirst + 1; i < ilast; i++)
{
dVc = current_cells[i].R * (current_x - current_cells[i].dif);
if ((dV_dcell && dVc * j_0e > 0 ||
(dV_dcell > 0 && (cell_data[i].potV + dVc) > (cell_data[count_cells + 1].potV)) ||
(dV_dcell < 0 && (cell_data[i].potV + dVc) < (cell_data[count_cells + 1].potV))))
//if (((dV_dcell && (dVc * j_0e > 0)) ||
// (dV_dcell > 0 && (cell_data[i].potV + dVc) > (cell_data[count_cells + 1].potV)) ||
// (dV_dcell < 0 && (cell_data[i].potV + dVc) < (cell_data[count_cells + 1].potV))))
if ((dV_dcell && (dVc * j_0e > 0)) ||
((dV_dcell > 0) && ((cell_data[i].potV + dVc) > cell_data[count_cells + 1].potV)) ||
((dV_dcell < 0) && ((cell_data[i].potV + dVc) < cell_data[count_cells + 1].potV)))
{
dVc = (cell_data[count_cells + 1].potV - cell_data[i].potV) / (count_cells + 1 - i);
}
cell_data[i + 1].potV = cell_data[i].potV + dVc;
}
if (!dV_dcell || fix_current)