diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index 8e375a4d..219d5a2b 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -3697,10 +3697,24 @@ factor(struct LOC_exec * LINK) n.UU.val = PhreeqcPtr->A0; break; 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; 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; case tokdh_av: n.UU.val = PhreeqcPtr->DH_Av; diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 1902f538..2c8a8541 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -1066,7 +1066,6 @@ public: LDBLE calc_vm_Cl(void); int multi_D(LDBLE DDt, int mobile_cell, 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); int fill_spec(int cell_no, int ref_cell); LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name); @@ -1693,7 +1692,7 @@ protected: int count_total_steps; 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, llnl_count_co2_coefs; diff --git a/phreeqcpp/Solution.cxx b/phreeqcpp/Solution.cxx index 3500128d..28e11fdf 100644 --- a/phreeqcpp/Solution.cxx +++ b/phreeqcpp/Solution.cxx @@ -1128,7 +1128,7 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble cxxNameDouble::iterator 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; } diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index 29be4855..6d824d5b 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -148,6 +148,7 @@ // #define MIN_RELATED_LOG_ACTIVITY -30 #endif #define REF_PRES_PASCAL 1.01325E5 /* Reference pressure: 1 atm */ +#define MAX_P_NONLLNL 1500.0 /* * Hash definitions */ diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 8b72bcb9..3055204f 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -237,6 +237,7 @@ initialize(void) /* Initialize llnl aqueous model parameters */ + a_llnl = b_llnl = 0.0; llnl_temp = (LDBLE *) PHRQ_malloc(sizeof(LDBLE)); if (llnl_temp == NULL) malloc_error(); diff --git a/phreeqcpp/model.cpp b/phreeqcpp/model.cpp index 1b8c2e46..e823bb4a 100644 --- a/phreeqcpp/model.cpp +++ b/phreeqcpp/model.cpp @@ -53,6 +53,11 @@ model(void) input_error++; 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) { @@ -563,7 +568,7 @@ gammas(LDBLE mu) */ int i, j; 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 muhalf, equiv; @@ -2935,7 +2940,7 @@ calc_gas_pressures(void) * Fixed-volume gas phase reacting with a solution * 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); 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); 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() + 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; //} - if (patm_x > 1500) - patm_x = 1500; + if (llnl_count_temp <= 0) + { + if (patm_x > MAX_P_NONLLNL) + patm_x = MAX_P_NONLLNL; + } } last_patm_x = patm_x; } diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp index 45446268..700809c1 100644 --- a/phreeqcpp/prep.cpp +++ b/phreeqcpp/prep.cpp @@ -5653,6 +5653,7 @@ calc_vm(LDBLE tc, LDBLE pa) * b4 = logk[vmi4], or * 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); for (int i = 0; i < count_s_x; i++) { diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp index 530090bf..49fc50b9 100644 --- a/phreeqcpp/print.cpp +++ b/phreeqcpp/print.cpp @@ -602,7 +602,7 @@ print_gas_phase(void) print_centered("Gas phase"); output_msg(sformatf("Total pressure: %5.2f atmospheres", (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"); else if (PR) output_msg(" (Peng-Robinson calculation)\n"); diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index e7001b20..918a86a3 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -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); } -/* ---------------------------------------------------------------------- */ -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:: 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]; 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; 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_j = A2; - if (sol_D[icell].tk_x == sol_D[jcell].tk_x) - b_j *= 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 * (f_free_j + g_j / ct[icell].visc2); + if (icell == count_cells && !stagnant) + 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 { - 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; + if (sol_D[icell].tk_x == sol_D[jcell].tk_x) + b_j *= sol_D[icell].spec[i].Dwt; + else + { + 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++; } @@ -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]; 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; 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; } } - b_i = A1; - b_j = A2 * sol_D[jcell].spec[j].Dwt; - if (sol_D[icell].tk_x == sol_D[jcell].tk_x) - b_i *= sol_D[jcell].spec[j].Dwt; + b_i = A1 * (f_free_i + g_i / ct[icell].visc1); + b_j = A2 * sol_D[jcell].spec[j].Dwt * (f_free_j + g_j / ct[icell].visc2); + if (icell == 0 && !stagnant) + 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 { - 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; + if (sol_D[icell].tk_x == sol_D[jcell].tk_x) + b_i *= sol_D[jcell].spec[j].Dwt; + else + { + 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++; } @@ -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; } } - b_i = A1 * sol_D[icell].spec[i].Dwt; - b_j = A2 * sol_D[jcell].spec[j].Dwt; - calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); + 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 * (f_free_j + g_j / ct[icell].visc2); + 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 //if (fabs(ddlm) > 1e-10) diff --git a/phreeqcpp/utilities.cpp b/phreeqcpp/utilities.cpp index 1a93255e..8ffc8c0b 100644 --- a/phreeqcpp/utilities.cpp +++ b/phreeqcpp/utilities.cpp @@ -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 + interpolation 0 - 300 oC, 0.006 - 1000 atm... */ + if (llnl_count_temp > 0) return OK; if (tc > 350.) { 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) + 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.) { tc = 350.;