max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 08:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Merge commit 'e1465e33228a3f867bd822405fb100614d5f6d0c'

Browse Source
This commit is contained in:
Darth Vader 2020-05-22 22:15:37 +00:00
commit 1c3edaa747
10 changed files with 107 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
phreeqcpp/PBasic.cpp
View File

@ -3697,10 +3697,24 @@ factor(struct LOC_exec * LINK)
n.UU.val = PhreeqcPtr->A0; n.UU.val = PhreeqcPtr->A0;
break; break;
case tokdh_a: case tokdh_a:
n.UU.val = PhreeqcPtr->DH_A; if (PhreeqcPtr->llnl_count_temp > 0)
{
n.UU.val = PhreeqcPtr->a_llnl;
}
else
{
n.UU.val = PhreeqcPtr->DH_A;
}
break; break;
case tokdh_b: case tokdh_b:
n.UU.val = PhreeqcPtr->DH_B; if (PhreeqcPtr->llnl_count_temp > 0)
{
n.UU.val = PhreeqcPtr->b_llnl;
}
else
{
n.UU.val = PhreeqcPtr->DH_B;
}
break; break;
case tokdh_av: case tokdh_av:
n.UU.val = PhreeqcPtr->DH_Av; n.UU.val = PhreeqcPtr->DH_Av;

3
phreeqcpp/Phreeqc.h
View File

@ -1066,7 +1066,6 @@ public:
LDBLE calc_vm_Cl(void); LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant); int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant); LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);
void calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant);
void diffuse_implicit(LDBLE DDt, int stagnant); void diffuse_implicit(LDBLE DDt, int stagnant);
int fill_spec(int cell_no, int ref_cell); int fill_spec(int cell_no, int ref_cell);
LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name); LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name);
@ -1693,7 +1692,7 @@ protected:
int count_total_steps; int count_total_steps;
int phast; int phast;
LDBLE *llnl_temp, *llnl_adh, *llnl_bdh, *llnl_bdot, *llnl_co2_coefs; LDBLE *llnl_temp, *llnl_adh, *llnl_bdh, *llnl_bdot, *llnl_co2_coefs, a_llnl, b_llnl;
int llnl_count_temp, llnl_count_adh, llnl_count_bdh, llnl_count_bdot, int llnl_count_temp, llnl_count_adh, llnl_count_bdh, llnl_count_bdot,
llnl_count_co2_coefs; llnl_count_co2_coefs;

2
phreeqcpp/Solution.cxx
View File

@ -1128,7 +1128,7 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble
cxxNameDouble::iterator it; cxxNameDouble::iterator it;
for (it = this->totals.begin(); it != this->totals.end(); it++) for (it = this->totals.begin(); it != this->totals.end(); it++)
{ {
if (it->second < 1e-18) if (it->second < 1e-25)
{ {
it->second = 0.0; it->second = 0.0;
} }

1
phreeqcpp/global_structures.h
View File

@ -148,6 +148,7 @@
// #define MIN_RELATED_LOG_ACTIVITY -30 // #define MIN_RELATED_LOG_ACTIVITY -30
#endif #endif
#define REF_PRES_PASCAL 1.01325E5 /* Reference pressure: 1 atm */ #define REF_PRES_PASCAL 1.01325E5 /* Reference pressure: 1 atm */
#define MAX_P_NONLLNL 1500.0
/* /*
* Hash definitions * Hash definitions
*/ */

1
phreeqcpp/mainsubs.cpp
View File

@ -237,6 +237,7 @@ initialize(void)
/* /*
Initialize llnl aqueous model parameters Initialize llnl aqueous model parameters
*/ */
a_llnl = b_llnl = 0.0;
llnl_temp = (LDBLE *) PHRQ_malloc(sizeof(LDBLE)); llnl_temp = (LDBLE *) PHRQ_malloc(sizeof(LDBLE));
if (llnl_temp == NULL) if (llnl_temp == NULL)
malloc_error(); malloc_error();

20
phreeqcpp/model.cpp
View File

@ -53,6 +53,11 @@ model(void)
input_error++; input_error++;
error_msg("Cannot use PITZER and SIT data blocks in same run (database + input file).", STOP); error_msg("Cannot use PITZER and SIT data blocks in same run (database + input file).", STOP);
} }
if ((pitzer_model == TRUE || sit_model == TRUE) && llnl_count_temp >0)
{
input_error++;
error_msg("Cannot use LLNL_AQUEOUS_MODEL_PARAMETERS with PITZER or SIT data blocks in same run (database + input file).", STOP);
}
if (pitzer_model == TRUE) if (pitzer_model == TRUE)
{ {
@ -563,7 +568,7 @@ gammas(LDBLE mu)
*/ */
int i, j; int i, j;
int ifirst, ilast; int ifirst, ilast;
LDBLE f, a_llnl, b_llnl, bdot_llnl, log_g_co2, dln_g_co2, c2_llnl; LDBLE f, bdot_llnl, log_g_co2, dln_g_co2, c2_llnl;
LDBLE c1, c2, a, b; LDBLE c1, c2, a, b;
LDBLE muhalf, equiv; LDBLE muhalf, equiv;
@ -2935,7 +2940,7 @@ calc_gas_pressures(void)
* Fixed-volume gas phase reacting with a solution * Fixed-volume gas phase reacting with a solution
* Change pressure used in logK to pressure of gas phase * Change pressure used in logK to pressure of gas phase
*/ */
if (gas_phase_ptr->Get_total_p() > 1500) if (gas_phase_ptr->Get_total_p() > MAX_P_NONLLNL && llnl_count_temp <= 0)
{ {
gas_phase_ptr->Set_total_moles(0); gas_phase_ptr->Set_total_moles(0);
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++) for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
@ -2945,12 +2950,12 @@ calc_gas_pressures(void)
struct phase *phase_ptr = phase_bsearch(gas_comp->Get_phase_name().c_str(), &j, FALSE); struct phase *phase_ptr = phase_bsearch(gas_comp->Get_phase_name().c_str(), &j, FALSE);
if (phase_ptr->in == TRUE) if (phase_ptr->in == TRUE)
{ {
phase_ptr->moles_x *= 1500.0 / gas_phase_ptr->Get_total_p(); phase_ptr->moles_x *= MAX_P_NONLLNL / gas_phase_ptr->Get_total_p();
gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() +
phase_ptr->moles_x); phase_ptr->moles_x);
} }
} }
gas_phase_ptr->Set_total_p(1500.0); gas_phase_ptr->Set_total_p(MAX_P_NONLLNL);
} }
} }
@ -3931,8 +3936,11 @@ reset(void)
//{ //{
// patm_x = ( 1 * patm_x + p_sat) / 2.0; // patm_x = ( 1 * patm_x + p_sat) / 2.0;
//} //}
if (patm_x > 1500) if (llnl_count_temp <= 0)
patm_x = 1500; {
if (patm_x > MAX_P_NONLLNL)
patm_x = MAX_P_NONLLNL;
}
} }
last_patm_x = patm_x; last_patm_x = patm_x;
} }

1
phreeqcpp/prep.cpp
View File

@ -5653,6 +5653,7 @@ calc_vm(LDBLE tc, LDBLE pa)
* b4 = logk[vmi4], or * b4 = logk[vmi4], or
* coef(tc) = millero[3] + millero[4] * tc + millero[5] * tc^2 * coef(tc) = millero[3] + millero[4] * tc + millero[5] * tc^2
*/ */
if (llnl_count_temp > 0) return OK;
LDBLE pb_s = 2600. + pa * 1.01325, TK_s = tc + 45.15, sqrt_mu = sqrt(mu_x); LDBLE pb_s = 2600. + pa * 1.01325, TK_s = tc + 45.15, sqrt_mu = sqrt(mu_x);
for (int i = 0; i < count_s_x; i++) for (int i = 0; i < count_s_x; i++)
{ {

2
phreeqcpp/print.cpp
View File

@ -602,7 +602,7 @@ print_gas_phase(void)
print_centered("Gas phase"); print_centered("Gas phase");
output_msg(sformatf("Total pressure: %5.2f atmospheres", output_msg(sformatf("Total pressure: %5.2f atmospheres",
(double) gas_phase_ptr->Get_total_p())); (double) gas_phase_ptr->Get_total_p()));
if (gas_phase_ptr->Get_total_p() >= 1500) if (gas_phase_ptr->Get_total_p() >= MAX_P_NONLLNL && llnl_count_temp <= 0)
output_msg(" WARNING: Program limit.\n"); output_msg(" WARNING: Program limit.\n");
else if (PR) else if (PR)
output_msg(" (Peng-Robinson calculation)\n"); output_msg(" (Peng-Robinson calculation)\n");

132
phreeqcpp/transport.cpp
View File

@ -3664,46 +3664,6 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant)
} }
return (OK); return (OK);
} }
/* ---------------------------------------------------------------------- */
void Phreeqc::
calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant)
/* ---------------------------------------------------------------------- */
{
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j));
// At filterends, concentrations of ions change step-wise to the DL.
// We take the harmonic mean for f_free, the average for the DL.
if (ct[icell].v_m[k].z)
{
if (!g_i && g_j)
{
ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) +
b_i * (1 - free_j) / 4 + b_j * g_j / 4;
}
else if (g_i && !g_j)
ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) +
b_j * (1 - free_i) / 4 + b_i * g_i / 4;
}
// for boundary cells...
if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2;
}
else
{
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j);
else if (icell == count_cells && jcell == count_cells + 1)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i);
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
return;
}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
LDBLE Phreeqc:: LDBLE Phreeqc::
find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
@ -4182,7 +4142,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else else
{ {
dum1 = it_sc->Get_mass_water() / t_aq2; dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1; dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_j += pow(dum2, ct[icell].v_m[k].z) * dum1; g_j += pow(dum2, ct[icell].v_m[k].z) * dum1;
} }
@ -4192,18 +4152,32 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
} }
} }
b_i = A1 * sol_D[icell].spec[i].Dwt; b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
b_j = A2; b_j = A2 * (f_free_j + g_j / ct[icell].visc2);
if (sol_D[icell].tk_x == sol_D[jcell].tk_x) if (icell == count_cells && !stagnant)
b_j *= sol_D[icell].spec[i].Dwt; ct[icell].v_m[k].b_ij = b_i;
else if (icell == all_cells - 1 && stagnant)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf *= 2 for this 'reservoir' cell in the input */
else else
{ {
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f; if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x); b_j *= sol_D[icell].spec[i].Dwt;
dum2 *= sol_D[jcell].viscos_f; else
b_j *= dum2; {
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
dum2 *= sol_D[jcell].viscos_f;
b_j *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == 0 && !stagnant)
ct[icell].v_m[k].b_ij = b_j;
else if (icell == 3 && stagnant && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for stagnant cell 3 in the input */
} }
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
k++; k++;
} }
@ -4275,7 +4249,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z]; g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else else
{ {
dum1 = it_sc->Get_mass_water() / t_aq1; dum1 = it_sc->Get_mass_water() / mass_water_bulk_x;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1; dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_i += pow(dum2, ct[icell].v_m[k].z) * dum1; g_i += pow(dum2, ct[icell].v_m[k].z) * dum1;
} }
@ -4292,18 +4266,31 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j *= sol_D[jcell].spec[j].erm_ddl; g_j *= sol_D[jcell].spec[j].erm_ddl;
} }
} }
b_i = A1; b_i = A1 * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * sol_D[jcell].spec[j].Dwt; b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
if (sol_D[icell].tk_x == sol_D[jcell].tk_x) if (icell == 0 && !stagnant)
b_i *= sol_D[jcell].spec[j].Dwt; ct[icell].v_m[k].b_ij = b_j;
else if (icell == 3 && stagnant && g_j && !g_i)
ct[icell].v_m[k].b_ij = b_j / 2; /* with the mixf *= 2 for 'reservoir' cell 3 in the input */
else else
{ {
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f; if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x); b_i *= sol_D[jcell].spec[j].Dwt;
dum2 *= sol_D[icell].viscos_f; else
b_i *= dum2; {
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
dum2 *= sol_D[icell].viscos_f;
b_i *= dum2;
}
ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
if (icell == count_cells && !stagnant)
ct[icell].v_m[k].b_ij = b_i;
else if (jcell == all_cells - 1 && stagnant && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2; /* with the mixf * 2 for this 'reservoir' cell in the input */
} }
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
k++; k++;
} }
@ -4386,9 +4373,28 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
g_j *= sol_D[jcell].spec[j].erm_ddl; g_j *= sol_D[jcell].spec[j].erm_ddl;
} }
} }
b_i = A1 * sol_D[icell].spec[i].Dwt; b_i = A1 * sol_D[icell].spec[i].Dwt * (f_free_i + g_i / ct[icell].visc1);
b_j = A2 * sol_D[jcell].spec[j].Dwt; b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2);
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); ct[icell].v_m[k].b_ij = b_i * b_j / (b_i + b_j);
// but for boundary cells...
if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i / 2;
}
else
{
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
ct[icell].v_m[k].b_ij = b_j;
else if (icell == count_cells && jcell == count_cells + 1)
ct[icell].v_m[k].b_ij = b_i;
}
if (ct[icell].v_m[k].z)
ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z;
//ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; // appt: this could give an incorrect large factor for implicit //ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; // appt: this could give an incorrect large factor for implicit
//if (fabs(ddlm) > 1e-10) //if (fabs(ddlm) > 1e-10)

2
phreeqcpp/utilities.cpp
View File

@ -167,6 +167,7 @@ calc_rho_0(LDBLE tc, LDBLE pa)
Wagner and Pruss, 2002, JPCRD 31, 387, eqn. 2.6, along the saturation pressure line + Wagner and Pruss, 2002, JPCRD 31, 387, eqn. 2.6, along the saturation pressure line +
interpolation 0 - 300 oC, 0.006 - 1000 atm... interpolation 0 - 300 oC, 0.006 - 1000 atm...
*/ */
if (llnl_count_temp > 0) return OK;
if (tc > 350.) if (tc > 350.)
{ {
if (need_temp_msg < 1) if (need_temp_msg < 1)
@ -226,6 +227,7 @@ calc_dielectrics(LDBLE tc, LDBLE pa)
and Fernandez et al., 1997, JPCRD 26, 1125, show its correctness) and Fernandez et al., 1997, JPCRD 26, 1125, show its correctness)
+ d(eps)/d(P), Debye-Hueckel A and B, and Av (for Av, see Pitzer et al., 1984, JPCRD 13, p. 4) + d(eps)/d(P), Debye-Hueckel A and B, and Av (for Av, see Pitzer et al., 1984, JPCRD 13, p. 4)
*/ */
if (llnl_count_temp > 0) return OK;
if (tc > 350.) if (tc > 350.)
{ {
tc = 350.; tc = 350.;