diff --git a/phreeqcpp/ExchComp.cxx b/phreeqcpp/ExchComp.cxx index dba46b3d..48d89dde 100644 --- a/phreeqcpp/ExchComp.cxx +++ b/phreeqcpp/ExchComp.cxx @@ -406,7 +406,7 @@ cxxExchComp::multiply(LDBLE extensive) { this->totals.multiply(extensive); this->charge_balance *= extensive; - this->phase_proportion *= extensive; + //this->phase_proportion *= extensive; } void diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index 219d5a2b..6f863c44 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -1612,6 +1612,9 @@ listtokens(FILE * f, tokenrec * l_buf) case tokdh_a: output_msg("DH_A"); // Debye-Hueckel A break; + case tokdebye_length: + output_msg("DEBYE_LENGTH"); // Debye-Hueckel length + break; case tokdh_b: output_msg("DH_B"); // Debye-Hueckel B break; @@ -3706,6 +3709,13 @@ factor(struct LOC_exec * LINK) n.UU.val = PhreeqcPtr->DH_A; } break; + case tokdebye_length: + { + double debye_length = (PhreeqcPtr->eps_r * EPSILON_ZERO * R_KJ_DEG_MOL * 1000.0 * PhreeqcPtr->tk_x) + / (2. * F_C_MOL * F_C_MOL * PhreeqcPtr->mu_x * 1000.); + n.UU.val = sqrt(debye_length); + break; + } case tokdh_b: if (PhreeqcPtr->llnl_count_temp > 0) { @@ -3993,7 +4003,7 @@ factor(struct LOC_exec * LINK) // call callback Basic function n.UU.val = (parse_all) ? 1 : PhreeqcPtr->basic_callback(x1, x2, str); - + PhreeqcPtr->PHRQ_free(str); } break; @@ -7377,6 +7387,7 @@ const std::map::value_type temp_tokens[] std::map::value_type("eps_r", PBasic::tokeps_r), std::map::value_type("vm", PBasic::tokvm), std::map::value_type("dh_a", PBasic::tokdh_a), + std::map::value_type("debye_length", PBasic::tokdebye_length), std::map::value_type("dh_b", PBasic::tokdh_b), std::map::value_type("dh_av", PBasic::tokdh_av), std::map::value_type("qbrn", PBasic::tokqbrn), diff --git a/phreeqcpp/PBasic.h b/phreeqcpp/PBasic.h index 76728e58..870e5850 100644 --- a/phreeqcpp/PBasic.h +++ b/phreeqcpp/PBasic.h @@ -317,6 +317,7 @@ public: tokvm, tokphase_vm, tokdh_a, + tokdebye_length, tokdh_b, tokdh_av, tokqbrn, diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 2c8a8541..b1089dfd 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -1033,15 +1033,19 @@ public: int tidy_logk(void); int tidy_exchange(void); int tidy_min_exchange(void); + int update_min_exchange(void); int tidy_kin_exchange(void); + int update_kin_exchange(void); int tidy_gas_phase(void); int tidy_inverse(void); int tidy_isotopes(void); int tidy_isotope_ratios(void); int tidy_isotope_alphas(void); int tidy_kin_surface(void); + int update_kin_surface(void); int tidy_master_isotope(void); int tidy_min_surface(void); + int update_min_surface(void); int tidy_phases(void); int tidy_pp_assemblage(void); int tidy_solutions(void); @@ -1066,6 +1070,7 @@ 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); diff --git a/phreeqcpp/class_main.cpp b/phreeqcpp/class_main.cpp index 2ee57fc6..0ee58bd1 100644 --- a/phreeqcpp/class_main.cpp +++ b/phreeqcpp/class_main.cpp @@ -109,7 +109,10 @@ main_method(int argc, char *argv[]) { return errors; } +#ifndef NO_UTF8_ENCODING #ifdef DOS + SetConsoleOutputCP(CP_UTF8); +#endif write_banner(); #endif @@ -199,7 +202,10 @@ main_method(int argc, char *argv[]) { return errors; } +#ifndef NO_UTF8_ENCODING #ifdef DOS + SetConsoleOutputCP(CP_UTF8); +#endif write_banner(); #endif @@ -271,9 +277,9 @@ write_banner(void) char buffer[80]; int len, indent; screen_msg( - " ÛßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÛ\n"); + " █▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀█\n"); screen_msg( - " º º\n"); + " â•‘ â•‘\n"); /* version */ #ifdef NPP @@ -282,21 +288,21 @@ write_banner(void) len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@"); #endif indent = (44 - len) / 2; - screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer, + screen_msg(sformatf("%14câ•‘%*c%s%*câ•‘\n", ' ', indent, ' ', buffer, 44 - indent - len, ' ')); screen_msg( - " º º\n"); + " â•‘ â•‘\n"); screen_msg( - " º A hydrogeochemical transport model º\n"); + " â•‘ A hydrogeochemical transport model â•‘\n"); screen_msg( - " º º\n"); + " â•‘ â•‘\n"); screen_msg( - " º by º\n"); + " â•‘ by â•‘\n"); screen_msg( - " º D.L. Parkhurst and C.A.J. Appelo º\n"); + " â•‘ D.L. Parkhurst and C.A.J. Appelo â•‘\n"); screen_msg( - " º º\n"); + " â•‘ â•‘\n"); /* date */ @@ -306,11 +312,11 @@ write_banner(void) len = sprintf(buffer, "%s", "@VER_DATE@"); #endif indent = (44 - len) / 2; - screen_msg(sformatf("%14cº%*c%s%*cº\n", ' ', indent, ' ', buffer, + screen_msg(sformatf("%14câ•‘%*c%s%*câ•‘\n", ' ', indent, ' ', buffer, 44 - indent - len, ' ')); screen_msg( - " ÛÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÛ\n\n"); + " █▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄█\n\n"); return 0; } @@ -485,7 +491,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, } local_database_file->close(); delete local_database_file; - + user_database = (char *) free_check_null(user_database); user_database = string_duplicate(token); screen_msg(sformatf("Database file: %s\n\n", token)); diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 3055204f..27e34981 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -2679,6 +2679,11 @@ save_init(int i) void Phreeqc::do_mixes(void) { + bool surf, exch, kin, min; + surf = (Rxn_surface_mix_map.size() > 0); + exch = (Rxn_exchange_mix_map.size() > 0); + kin = (Rxn_kinetics_mix_map.size() > 0); + min = (Rxn_pp_assemblage_mix_map.size() > 0); Utilities::Rxn_mix(Rxn_solution_mix_map, Rxn_solution_map, this); Utilities::Rxn_mix(Rxn_exchange_mix_map, Rxn_exchange_map, this); Utilities::Rxn_mix(Rxn_gas_phase_mix_map, Rxn_gas_phase_map, this); @@ -2686,4 +2691,8 @@ Phreeqc::do_mixes(void) Utilities::Rxn_mix(Rxn_pp_assemblage_mix_map, Rxn_pp_assemblage_map, this); Utilities::Rxn_mix(Rxn_ss_assemblage_mix_map, Rxn_ss_assemblage_map, this); Utilities::Rxn_mix(Rxn_surface_mix_map, Rxn_surface_map, this); + if (exch || kin) update_kin_exchange(); + if (exch || min) update_min_exchange(); + if (surf || min) update_min_surface(); + if (surf || kin) update_kin_surface(); } diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp index 49fc50b9..1b8ab957 100644 --- a/phreeqcpp/print.cpp +++ b/phreeqcpp/print.cpp @@ -1479,11 +1479,12 @@ print_species(void) { output_msg(sformatf("%50s%10s%10s%10s\n", "Log", "Log", "Log", "mole V")); } - output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%10s\n\n", "Species", #ifdef NO_UTF8_ENCODING + output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%10s\n\n", "Species", "Molality", "Activity", "Molality", "Activity", "Gamma", "cm3/mol")); #else - "Molality", "Activity", "Molality", "Activity", "Gamma", "cm³/mol")); + output_msg(sformatf(" %-13s%12s%12s%10s%10s%10s%11s\n\n", "Species", + "Molality", "Activity", "Molality", "Activity", "Gamma", "cm³/mol")); #endif /* * Print list of species @@ -1651,7 +1652,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING output_msg(sformatf("\t%11.3e sigma, C/m2\n", #else - output_msg(sformatf("\t%11.3e sigma, C/m²\n", + output_msg(sformatf("\t%11.3e sigma, C/m²\n", #endif (double) (charge * F_C_MOL / (charge_ptr->Get_specific_area() * @@ -1662,7 +1663,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING output_msg(sformatf("\tundefined sigma, C/m2\n")); #else - output_msg(sformatf("\tundefined sigma, C/m²\n")); + output_msg(sformatf("\tundefined sigma, C/m²\n")); #endif } if (use.Get_surface_ptr()->Get_type() == cxxSurface::CCM) @@ -1684,7 +1685,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING "\t%11.3e specific area, m2/mol %s\n", #else - "\t%11.3e specific area, m²/mol %s\n", + "\t%11.3e specific area, m²/mol %s\n", #endif (double) charge_ptr->Get_specific_area(), comp_ptr->Get_phase_name().c_str())); @@ -1692,7 +1693,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING "\t%11.3e m2 for %11.3e moles of %s\n\n", #else - "\t%11.3e m² for %11.3e moles of %s\n\n", + "\t%11.3e m² for %11.3e moles of %s\n\n", #endif (double) (charge_ptr->Get_grams() * charge_ptr->Get_specific_area()), @@ -1705,7 +1706,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING "\t%11.3e specific area, m2/mol %s\n", #else - "\t%11.3e specific area, m²/mol %s\n", + "\t%11.3e specific area, m²/mol %s\n", #endif (double) charge_ptr->Get_specific_area(), comp_ptr->Get_rate_name().c_str())); @@ -1713,7 +1714,7 @@ print_surface(void) #ifdef NO_UTF8_ENCODING "\t%11.3e m2 for %11.3e moles of %s\n\n", #else - "\t%11.3e m² for %11.3e moles of %s\n\n", + "\t%11.3e m² for %11.3e moles of %s\n\n", #endif (double) (charge_ptr->Get_grams() * charge_ptr->Get_specific_area()), @@ -1726,13 +1727,13 @@ print_surface(void) #ifdef NO_UTF8_ENCODING "\t%11.3e specific area, m2/g\n", #else - "\t%11.3e specific area, m²/g\n", + "\t%11.3e specific area, m²/g\n", #endif (double) charge_ptr->Get_specific_area())); #ifdef NO_UTF8_ENCODING output_msg(sformatf("\t%11.3e m2 for %11.3e g\n\n", #else - output_msg(sformatf("\t%11.3e m² for %11.3e g\n\n", + output_msg(sformatf("\t%11.3e m² for %11.3e g\n\n", #endif (double) (charge_ptr->Get_specific_area() * charge_ptr->Get_grams()), @@ -1948,28 +1949,28 @@ print_surface_cd_music(void) #ifdef NO_UTF8_ENCODING "\t%11.3e sigma, plane 0, C/m2\n", #else - "\t%11.3e sigma, plane 0, C/m²\n", + "\t%11.3e sigma, plane 0, C/m²\n", #endif (double) charge_ptr->Get_sigma0())); output_msg(sformatf( #ifdef NO_UTF8_ENCODING "\t%11.3e sigma, plane 1, C/m2\n", #else - "\t%11.3e sigma, plane 1, C/m²\n", + "\t%11.3e sigma, plane 1, C/m²\n", #endif (double) charge_ptr->Get_sigma1())); output_msg(sformatf( #ifdef NO_UTF8_ENCODING "\t%11.3e sigma, plane 2, C/m2\n", #else - "\t%11.3e sigma, plane 2, C/m²\n", + "\t%11.3e sigma, plane 2, C/m²\n", #endif (double) charge_ptr->Get_sigma2())); output_msg(sformatf( #ifdef NO_UTF8_ENCODING "\t%11.3e sigma, diffuse layer, C/m2\n\n", #else - "\t%11.3e sigma, diffuse layer, C/m²\n\n", + "\t%11.3e sigma, diffuse layer, C/m²\n\n", #endif (double) charge_ptr->Get_sigmaddl())); } @@ -1978,7 +1979,7 @@ print_surface_cd_music(void) #ifdef NO_UTF8_ENCODING output_msg(sformatf("\tundefined sigma, C/m2\n")); #else - output_msg(sformatf("\tundefined sigma, C/m²\n")); + output_msg(sformatf("\tundefined sigma, C/m²\n")); #endif } output_msg(sformatf("\t%11.3e psi, plane 0, V\n", @@ -2232,11 +2233,12 @@ print_totals(void) if (SC > 0) { //output_msg(sformatf("%36s%i%7s%i\n", - output_msg(sformatf("%35s%3.0f%7s%i\n", #ifdef NO_UTF8_ENCODING + output_msg(sformatf("%35s%3.0f%7s%i\n", "Specific Conductance (uS/cm, ", tc_x, "oC) = ", (int) SC)); #else - "Specific Conductance (µS/cm, ", tc_x, "°C) = ", (int) SC)); + output_msg(sformatf("%36s%3.0f%7s%i\n", + "Specific Conductance (µS/cm, ", tc_x, "°C) = ", (int) SC)); #endif } /* VP: Density Start */ @@ -2246,7 +2248,7 @@ print_totals(void) #ifdef NO_UTF8_ENCODING output_msg(sformatf("%45s%9.5f", "Density (g/cm3) = ", #else - output_msg(sformatf("%45s%9.5f", "Density (g/cm³) = ", + output_msg(sformatf("%46s%9.5f", "Density (g/cm³) = ", #endif (double) dens)); if (state == INITIAL_SOLUTION && use.Get_solution_ptr()->Get_initial_data()->Get_calc_density()) @@ -2266,11 +2268,12 @@ print_totals(void) (double) viscos)); if (tc_x > 200 && !pure_water) { - output_msg(sformatf("%18s\n", #ifdef NO_UTF8_ENCODING + output_msg(sformatf("%18s\n", " (solute contributions limited to 200 oC)")); #else - " (solute contributions limited to 200 °C)")); + output_msg(sformatf("%19s\n", + " (solute contributions limited to 200 °C)")); #endif } else output_msg(sformatf("\n")); @@ -2300,7 +2303,7 @@ print_totals(void) #ifdef NO_UTF8_ENCODING output_msg(sformatf("%45s%6.2f\n", "Temperature (oC) = ", #else - output_msg(sformatf("%45s%6.2f\n", "Temperature (°C) = ", + output_msg(sformatf("%46s%6.2f\n", "Temperature (°C) = ", #endif (double) tc_x)); diff --git a/phreeqcpp/readtr.cpp b/phreeqcpp/readtr.cpp index 86b74682..e45f7a37 100644 --- a/phreeqcpp/readtr.cpp +++ b/phreeqcpp/readtr.cpp @@ -915,9 +915,9 @@ read_transport(void) warning_msg(error_string); for (i = count_por; i < all_cells - st; i++) { - if (i == max_cells) - continue; - assert((i+1) < all_cells); + //if (i == max_cells) + // continue; + //assert((i+1) < all_cells); if ((i+1) < all_cells) { cell_data[i + 1].por = pors[count_por - 1]; diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp index 09502d18..95d12e21 100644 --- a/phreeqcpp/tidy.cpp +++ b/phreeqcpp/tidy.cpp @@ -292,7 +292,45 @@ tidy_model(void) { tidy_solutions(); } - +/* +* need to update exchange and surface related in case anything has changed +*/ + if (keycount[Keywords::KEY_KINETICS] > 0 || + keycount[Keywords::KEY_KINETICS_RAW] > 0 || + keycount[Keywords::KEY_KINETICS_MODIFY] || + keycount[Keywords::KEY_EXCHANGE] > 0 || + keycount[Keywords::KEY_EXCHANGE_RAW] > 0 || + keycount[Keywords::KEY_EXCHANGE_MODIFY]) + { + update_kin_exchange(); + } + if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 || + keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 || + keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] || + keycount[Keywords::KEY_EXCHANGE] > 0 || + keycount[Keywords::KEY_EXCHANGE_RAW] > 0 || + keycount[Keywords::KEY_EXCHANGE_MODIFY]) + { + update_min_exchange(); + } + if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 || + keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 || + keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] || + keycount[Keywords::KEY_SURFACE] > 0 || + keycount[Keywords::KEY_SURFACE_RAW] > 0 || + keycount[Keywords::KEY_SURFACE_MODIFY] > 0) + { + update_min_surface(); + } + if (keycount[Keywords::KEY_KINETICS] > 0 || + keycount[Keywords::KEY_KINETICS_RAW] > 0 || + keycount[Keywords::KEY_KINETICS_MODIFY] > 0 || + keycount[Keywords::KEY_SURFACE] > 0 || + keycount[Keywords::KEY_SURFACE_RAW] > 0 || + keycount[Keywords::KEY_SURFACE_MODIFY] > 0) + { + update_kin_surface(); + } /* if (new_model || new_exchange || new_pp_assemblage || new_surface || new_gas_phase || new_kinetics) reset_last_model(); */ if (new_model) { @@ -3575,10 +3613,10 @@ tidy_kin_exchange(void) assert(false); } cxxExchange * exchange_ptr = &(it->second); - //if (!exchange_ptr->Get_new_def()) - // continue; - //if (exchange_ptr->Get_n_user() < 0) - // continue; + if (!exchange_ptr->Get_new_def()) + continue; + if (exchange_ptr->Get_n_user() < 0) + continue; // check elements for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++) { @@ -3669,6 +3707,140 @@ tidy_kin_exchange(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +update_kin_exchange(void) +/* ---------------------------------------------------------------------- */ +/* + * If exchanger is related to mineral, exchanger amount is + * set in proportion. Exchange needs to be updated if the + * amount of kinetic reaction has changed. Corner case of + * zero moles. + */ +{ + cxxKinetics* kinetics_ptr; + char* ptr; + LDBLE conc; + + std::map::iterator it = Rxn_exchange_map.begin(); + for ( ; it != Rxn_exchange_map.end(); it++) + { + cxxExchange* exchange_ptr = &(it->second); + if (exchange_ptr->Get_n_user() < 0) continue; + // check elements + for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++) + { + cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j]; + if (comp_ref.Get_rate_name().size() == 0) continue; + double comp_moles = 0.0; + /* First find exchange master species */ + cxxNameDouble nd = comp_ref.Get_totals(); + cxxNameDouble::iterator kit = nd.begin(); + bool found_exchange = false; + for (; kit != nd.end(); kit++) + { + /* Find master species */ + struct element* elt_ptr = element_store(kit->first.c_str()); + if (elt_ptr == NULL || elt_ptr->master == NULL) + { + input_error++; + error_string = sformatf("Master species not in database " + "for %s, skipping element.", + kit->first.c_str()); + error_msg(error_string, CONTINUE); + continue; + } + if (elt_ptr->master->type == EX) + { + comp_moles = kit->second; + found_exchange = true; + } + } + //if (!found_exchange) + //{ + // input_error++; + // error_string = sformatf( + // "Exchange formula does not contain an exchange master species, %s", + // comp_ref.Get_formula().c_str()); + // error_msg(error_string, CONTINUE); + // continue; + //} + + /* Now find associated kinetic reaction ... */ + if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, exchange_ptr->Get_n_user())) == NULL) + { + input_error++; + error_string = sformatf( + "Kinetics %d must be defined to use exchange related to kinetic reaction, %s", + exchange_ptr->Get_n_user(), comp_ref.Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + size_t k; + for (k = 0; k < kinetics_ptr->Get_kinetics_comps().size(); k++) + { + if (strcmp_nocase + (comp_ref.Get_rate_name().c_str(), + kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str()) == 0) + { + break; + } + } + if (k == kinetics_ptr->Get_kinetics_comps().size()) + { + input_error++; + error_string = sformatf( + "Kinetic reaction, %s, related to exchanger, %s, not found in KINETICS %d", + comp_ref.Get_rate_name().c_str(), comp_ref.Get_formula().c_str(), exchange_ptr->Get_n_user()); + error_msg(error_string, CONTINUE); + continue; + } + + /* use database name for phase */ + comp_ref.Set_rate_name(kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str()); + + /* make exchanger concentration proportional to mineral ... */ + conc = kinetics_ptr->Get_kinetics_comps()[k].Get_m() * comp_ref.Get_phase_proportion(); + if (found_exchange && comp_moles > 0.0) + { + /* parse formula */ + count_elts = 0; + paren_count = 0; + { + char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str()); + ptr = temp_formula; + get_elts_in_species(&ptr, 1.0); + free_check_null(temp_formula); + } + cxxNameDouble nd_formula = elt_list_NameDouble(); + double comp_coef = 0; + for (kit = nd_formula.begin(); kit != nd_formula.end(); kit++) + { + /* Find master species */ + struct element* elt_ptr = element_store(kit->first.c_str()); + if (elt_ptr->master->type == EX) + { + comp_coef = kit->second; + } + } + comp_ref.multiply(comp_coef * conc / comp_moles); + } + else /* need to generate totals from scratch */ + { + count_elts = 0; + paren_count = 0; + { + char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str()); + ptr = temp_formula; + get_elts_in_species(&ptr, conc); + free_check_null(temp_formula); + } + comp_ref.Set_totals(elt_list_NameDouble()); + } + } + } + return (OK); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: tidy_min_exchange(void) /* ---------------------------------------------------------------------- */ /* @@ -3694,10 +3866,10 @@ tidy_min_exchange(void) assert(false); } cxxExchange * exchange_ptr = &(it->second); - //if (!exchange_ptr->Get_new_def()) - // continue; - //if (exchange_ptr->Get_n_user() < 0) - // continue; + if (!exchange_ptr->Get_new_def()) + continue; + if (exchange_ptr->Get_n_user() < 0) + continue; n = exchange_ptr->Get_n_user(); // check elements for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++) @@ -3832,6 +4004,187 @@ tidy_min_exchange(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +update_min_exchange(void) +/* ---------------------------------------------------------------------- */ +/* + * If exchanger is related to mineral, exchanger amount is + * set in proportion. Need to check in case exchange or min + * are modified. + */ +{ + int n, jj; + char* ptr; + LDBLE conc; + + std::map::iterator it = Rxn_exchange_map.begin(); + for ( ; it != Rxn_exchange_map.end(); it++) + { + cxxExchange* exchange_ptr = &(it->second); + if (exchange_ptr->Get_n_user() < 0) continue; + n = exchange_ptr->Get_n_user(); + // check elements + for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++) + { + double comp_moles = 0.0; + cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j]; + if (comp_ref.Get_phase_name().size() == 0) continue; + /* First find exchange master species */ + cxxNameDouble nd = comp_ref.Get_totals(); + cxxNameDouble::iterator kit = nd.begin(); + bool found_exchange = false; + for (; kit != nd.end(); kit++) + { + /* Find master species */ + struct element* elt_ptr = element_store(kit->first.c_str()); + if (elt_ptr == NULL || elt_ptr->master == NULL) + { + input_error++; + error_string = sformatf("Master species not in database " + "for %s, skipping element.", + kit->first.c_str()); + error_msg(error_string, CONTINUE); + continue; + } + if (elt_ptr->master->type == EX) + { + comp_moles = kit->second; + found_exchange = true; + } + } + //if (!found_exchange) + //{ + // input_error++; + // error_string = sformatf( + // "Exchange formula does not contain an exchange master species, %s", + // comp_ref.Get_formula().c_str()); + // error_msg(error_string, CONTINUE); + // continue; + //} + + cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n); + /* Now find the mineral on which exchanger depends... */ + if (pp_assemblage_ptr == NULL) + { + input_error++; + error_string = sformatf( + "Equilibrium_phases %d must be defined to use exchange related to mineral phase, %s", + n, comp_ref.Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + std::map::iterator jit; + jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin(); + for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++) + { + if (strcmp_nocase(comp_ref.Get_phase_name().c_str(), jit->first.c_str()) == 0) + { + break; + } + } + if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end()) + { + input_error++; + error_string = sformatf( + "Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d", + comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n); + error_msg(error_string, CONTINUE); + continue; + } + /* use database name for phase */ + comp_ref.Set_phase_name(jit->first.c_str()); + + /* make exchanger concentration proportional to mineral ... */ + conc = jit->second.Get_moles() * comp_ref.Get_phase_proportion(); + if (found_exchange && comp_moles > 0.0) + { + /* parse formula */ + count_elts = 0; + paren_count = 0; + { + char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str()); + ptr = temp_formula; + get_elts_in_species(&ptr, 1.0); + free_check_null(temp_formula); + } + cxxNameDouble nd_formula = elt_list_NameDouble(); + double comp_coef = 0; + for (kit = nd_formula.begin(); kit != nd_formula.end(); kit++) + { + /* Find master species */ + struct element* elt_ptr = element_store(kit->first.c_str()); + if (elt_ptr->master->type == EX) + { + comp_coef = kit->second; + } + } + comp_ref.multiply(comp_coef * conc / comp_moles); + } + else /* comp_moles is zero, need to redefine totals from scratch */ + { + count_elts = 0; + paren_count = 0; + { + char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str()); + ptr = temp_formula; + get_elts_in_species(&ptr, conc); + free_check_null(temp_formula); + } + comp_ref.Set_totals(elt_list_NameDouble()); + /* + * make sure exchange elements are in phase + */ + count_elts = 0; + paren_count = 0; + { + char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str()); + ptr = temp_formula; + get_elts_in_species(&ptr, -comp_ref.Get_phase_proportion()); + free_check_null(temp_formula); + } + int l; + struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE); + if (phase_ptr != NULL) + { + char* temp_formula = string_duplicate(phase_ptr->formula); + ptr = temp_formula; + get_elts_in_species(&ptr, 1.0); + free_check_null(temp_formula); + } + else + { + input_error++; + error_string = sformatf( + "Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d", + comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n); + error_msg(error_string, CONTINUE); + continue; + } + qsort(elt_list, (size_t)count_elts, + (size_t)sizeof(struct elt_list), elt_list_compare); + elt_list_combine(); + for (jj = 0; jj < count_elts; jj++) + { + if (elt_list[jj].elt->primary->s->type != EX + && elt_list[jj].coef < 0) + { + input_error++; + error_string = sformatf( + "Stoichiometry of exchanger, %s * %g mol sites/mol phase,\n\tmust be a subset of the related phase %s, %s.", + comp_ref.Get_formula().c_str(), + (double)comp_ref.Get_phase_proportion(), + phase_ptr->name, + phase_ptr->formula); + error_msg(error_string, CONTINUE); + break; + } + } + } + } + } + return (OK); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: tidy_min_surface(void) /* ---------------------------------------------------------------------- */ /* @@ -3853,7 +4206,6 @@ tidy_min_surface(void) assert(false); } cxxSurface *surface_ptr = &(kit->second); - if (!surface_ptr->Get_new_def()) continue; if (!surface_ptr->Get_new_def()) continue; if (surface_ptr->Get_n_user() < 0) @@ -4084,6 +4436,151 @@ tidy_min_surface(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +update_min_surface(void) +/* ---------------------------------------------------------------------- */ +/* + * If surface is related to mineral, surface amount is + * set in proportion + */ +{ + std::map::iterator kit; + for (kit = Rxn_surface_map.begin(); kit != Rxn_surface_map.end(); kit++) + { + cxxSurface* surface_ptr = &(kit->second); + if (surface_ptr->Get_n_user() < 0) continue; + for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++) + { + double comp_moles = 0.0; + cxxSurfaceComp* surface_comp_ptr = &(surface_ptr->Get_surface_comps()[j]); + if (surface_comp_ptr->Get_phase_name().size() == 0) continue; + cxxSurfaceCharge* surface_charge_ptr = NULL; + if (surface_ptr->Get_type() != cxxSurface::NO_EDL) + { + surface_charge_ptr = surface_ptr->Find_charge(surface_comp_ptr->Get_charge_name()); + if (surface_charge_ptr == NULL) + { + input_error++; + error_string = sformatf("Data structure for surface charge not found " + "for %s ", + surface_comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + } + int n = surface_ptr->Get_n_user(); + + /* First find surface master species */ + cxxNameDouble::iterator it; + for (it = surface_comp_ptr->Get_totals().begin(); it != surface_comp_ptr->Get_totals().end(); it++) + { + /* Find master species */ + struct element* elt_ptr = element_store(it->first.c_str()); + struct master* master_ptr = elt_ptr->master; + if (master_ptr == NULL) + { + input_error++; + error_string = sformatf("Master species not in database " + "for %s, skipping element.", + elt_ptr->name); + error_msg(error_string, CONTINUE); + continue; + } + if (master_ptr->type != SURF) continue; + comp_moles = it->second; + surface_comp_ptr->Set_master_element(elt_ptr->name); + break; + } + //if (surface_comp_ptr->Get_master_element().size() == 0) + //{ + // input_error++; + // error_string = sformatf( + // "Surface formula does not contain a surface master species, %s", + // surface_comp_ptr->Get_formula().c_str()); + // error_msg(error_string, CONTINUE); + // continue; + //} + + /* Now find the mineral on which surface depends... */ + cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n); + if (pp_assemblage_ptr == NULL) + { + input_error++; + error_string = sformatf( + "Equilibrium_phases %d must be defined to use surface related to mineral phase, %s", + n, surface_comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + std::map::iterator jit; + jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin(); + for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++) + { + if (strcmp_nocase(surface_comp_ptr->Get_phase_name().c_str(), + jit->first.c_str()) == 0) + { + break; + } + } + if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end()) + { + input_error++; + error_string = sformatf( + "Mineral, %s, related to surface, %s, not found in Equilibrium_Phases %d", + surface_comp_ptr->Get_phase_name().c_str(), surface_comp_ptr->Get_formula().c_str(), n); + error_msg(error_string, CONTINUE); + continue; + } + int l; + struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE); + if (phase_ptr == NULL) + { + input_error++; + error_string = sformatf( + "Mineral, %s, related to surface, %s, not found in database.", + jit->first.c_str(), surface_comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + /* use database name for phase */ + surface_comp_ptr->Set_phase_name(phase_ptr->name); + /* make surface concentration proportional to mineral ... */ + LDBLE conc = jit->second.Get_moles() * surface_comp_ptr->Get_phase_proportion(); + double grams = 0.0; + if (surface_charge_ptr != NULL) + { + grams = surface_charge_ptr->Get_grams(); + } + if (comp_moles > 0.0) + { + surface_comp_ptr->multiply(conc / comp_moles); + } + else /* need to generate from scratch */ + { + char* temp_formula = string_duplicate(surface_comp_ptr->Get_formula().c_str()); + char* ptr = temp_formula; + count_elts = 0; + paren_count = 0; + get_elts_in_species(&ptr, conc); + free_check_null(temp_formula); + + cxxNameDouble nd = elt_list_NameDouble(); + surface_comp_ptr->Set_totals(nd); + } + if (grams > 0.0) + { + surface_charge_ptr->multiply(jit->second.Get_moles() / grams); + } + else if (surface_charge_ptr != NULL) /* need to generate from scratch */ + { + surface_charge_ptr->Set_grams(jit->second.Get_moles()); + surface_charge_ptr->Set_charge_balance(0.0); + } + } + } + return (OK); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: tidy_kin_surface(void) /* ---------------------------------------------------------------------- */ /* @@ -4367,6 +4864,144 @@ tidy_kin_surface(void) } /* ---------------------------------------------------------------------- */ int Phreeqc:: +update_kin_surface(void) +/* ---------------------------------------------------------------------- */ +/* + * If surface is related to mineral, surface amount is + * set in proportion. Need to update surface if + * moles of kinetic reaction changes + */ +{ + cxxKinetics* kinetics_ptr; + + std::map::iterator it; + for (it = Rxn_surface_map.begin(); it != Rxn_surface_map.end(); it++) + { + cxxSurface* surface_ptr = &(it->second); + if (surface_ptr->Get_n_user() < 0) continue; + int n = surface_ptr->Get_n_user(); + for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++) + { + double comp_moles = 0.0; + cxxSurfaceComp* comp_ptr = &(surface_ptr->Get_surface_comps()[j]); + if (comp_ptr->Get_rate_name().size() == 0) continue; + comp_ptr->Set_master_element(""); + + /* First find surface master species */ + int k; + cxxNameDouble::iterator kit; + for (kit = comp_ptr->Get_totals().begin(); kit != comp_ptr->Get_totals().end(); kit++) + { + /* Find master species */ + struct element* elt_ptr = element_store(kit->first.c_str()); + struct master* master_ptr = elt_ptr->master; + if (master_ptr == NULL) + { + input_error++; + error_string = sformatf("Master species not in database " + "for %s, skipping element.", + elt_ptr->name); + error_msg(error_string, CONTINUE); + continue; + } + if (master_ptr->type != SURF) continue; + comp_ptr->Set_master_element(elt_ptr->name); + comp_moles = kit->second; + break; + } + if (comp_ptr->Get_master_element().size() == 0) + { + input_error++; + error_string = sformatf( + "Surface formula does not contain a surface master species, %s", + comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + + /* Now find the kinetic reaction on which surface depends... */ + if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n)) == NULL) + { + input_error++; + error_string = sformatf( + "Kinetics %d must be defined to use surface related to kinetic reaction, %s", + n, comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + for (k = 0; k < (int)kinetics_ptr->Get_kinetics_comps().size(); k++) + { + cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]); + if (strcmp_nocase + (comp_ptr->Get_rate_name().c_str(), + kin_comp_ptr->Get_rate_name().c_str()) == 0) + { + break; + } + } + if (k == (int)kinetics_ptr->Get_kinetics_comps().size()) + { + input_error++; + error_string = sformatf( + "Kinetic reaction, %s, related to surface, %s, not found in Kinetics %d", + comp_ptr->Get_rate_name().c_str(), comp_ptr->Get_formula().c_str(), n); + error_msg(error_string, CONTINUE); + continue; + } + + cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]); + /* use database name for rate */ + comp_ptr->Set_rate_name(kin_comp_ptr->Get_rate_name().c_str()); + cxxSurfaceCharge* charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name()); + if (surface_ptr->Get_type() != cxxSurface::NO_EDL) + { + charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name()); + if (charge_ptr == NULL) + { + input_error++; + error_string = sformatf("Data structure for surface charge not found " + "for %s ", + comp_ptr->Get_formula().c_str()); + error_msg(error_string, CONTINUE); + continue; + } + } + /* make surface concentration proportional to mineral ... */ + LDBLE conc = kin_comp_ptr->Get_m() * comp_ptr->Get_phase_proportion(); + double grams = 0.0; + if (charge_ptr != NULL) charge_ptr->Get_grams(); + if (comp_moles > 0.0) + { + comp_ptr->multiply(conc / comp_moles); + } + else /* need to generate from scratch */ + { + char* temp_formula = string_duplicate(comp_ptr->Get_formula().c_str()); + char* ptr = temp_formula; + count_elts = 0; + paren_count = 0; + get_elts_in_species(&ptr, conc); + free_check_null(temp_formula); + + cxxNameDouble nd = elt_list_NameDouble(); + comp_ptr->Set_totals(nd); + } + + if (grams > 0.0) + { + charge_ptr->multiply(kin_comp_ptr->Get_m() / grams); + } + else if (charge_ptr != NULL) /* need to generate from scratch */ + { + charge_ptr->Set_grams(kin_comp_ptr->Get_m()); + charge_ptr->Set_charge_balance(0.0); + } + } + } + return (OK); +} +/* ---------------------------------------------------------------------- */ +int Phreeqc:: ss_prep(LDBLE t, cxxSS *ss_ptr, int print) /* ---------------------------------------------------------------------- */ { diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 918a86a3..8684ba62 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -1,4 +1,4 @@ -#include "Utils.h" +#include "Utils.h" #include "Phreeqc.h" #include "phqalloc.h" #include "Exchange.h" @@ -585,7 +585,6 @@ transport(void) if (overall_iterations > max_iter) max_iter = overall_iterations; cell_no = i; - mixrun = j; if (multi_Dflag) sprintf(token, "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", @@ -746,7 +745,6 @@ transport(void) if (i == first_c && count_cells > 1) kin_time /= 2; cell_no = i; - mixrun = 0; if (multi_Dflag) sprintf(token, "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", @@ -794,6 +792,7 @@ transport(void) j = 1; for (; j <= nmix; j++) // loop on j { + mixrun = j; if (multi_Dflag && j == nmix && (transport_step % print_modulus == 0)) { sprintf(token, @@ -850,7 +849,6 @@ transport(void) if (overall_iterations > max_iter) max_iter = overall_iterations; cell_no = i; - mixrun = j; if (multi_Dflag) sprintf(token, "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", @@ -2827,7 +2825,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) } current_A = current_x / DDt * F_C_MOL; - for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++) + for (i = ifirst; i <= ilast + stagnant + ((bcon_last == 2 || (dV_dcell && stagnant)) ? 1 : 0); i++) { if (i <= ilast + 1) { @@ -2879,10 +2877,12 @@ diffuse_implicit(LDBLE DDt, int stagnant) for (icell = if1; icell != il1; icell += incr) { min_mol = min_dif_M * ct[icell].kgw; + if (min_mol < 1e-13) + min_mol = 1e-13; dum1 = dum2 = 0; sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell); sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1); - if (stagnant) + if (stagnant && mixf_stag[icell][cp]) { i1 = (icell == 0 ? c1 + 1 : icell == c1 ? cc1 : icell + c1); sptr_stag = Utilities::Rxn_find(Rxn_solution_map, i1); @@ -2890,6 +2890,13 @@ diffuse_implicit(LDBLE DDt, int stagnant) else sptr_stag = NULL; + //if (!cp) + //{ + // ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0; + // if (sptr_stag) + // ct[i1].J_ij_sum = 0.0; + //} + if (!strcmp(ct[icell].m_s[cp].name, "H")) { dummy = ct[icell].m_s[cp].tot1; @@ -2933,7 +2940,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) } sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy); } - //if (cp == count_m_s - 1) // transport the charge imbalance (but doesn't improve the change of H2 or O2). + //if (cp == count_m_s - 1) // transport the charge imbalance //{ // sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum); // sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum); @@ -2942,180 +2949,166 @@ diffuse_implicit(LDBLE DDt, int stagnant) //} continue; } - // see if icell - incr has negative moles, subtract it, so that it is canceled... - if (icell > 0 && icell <= ilast && - (it1 = neg_moles.find(icell - incr)) != neg_moles.end() - && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) - { - ct[icell].m_s[cp].tot1 += it2->second; - neg_moles.erase(it1); - els.erase(it2); - neg_moles.insert(std::make_pair(icell - incr, els)); - } dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name]; - if (!dum1) + dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name]; + if (sptr_stag) + dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name]; + + // check for negative moles, add moles from other redox states and the donnan layer when necessary and available... + if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol && + (dV_dcell || (icell > 0 && icell <= ilast))) { dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name); - if (dum1) - sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1; - } - dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name]; - if (!dum2) - { - dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name); - if (dum2) - sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2; - } - if (sptr_stag) - { - dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name]; - if (!dum_stag) - { - dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name); - if (dum_stag) - sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag; - } - } - - // check for negative moles, add moles from the donnan layer when necessary and available... - if (ct[icell].m_s[cp].tot1 > dum1 && - (dV_dcell || (icell > 0 && icell <= ilast))) - { - if (!ct[icell].dl_s) - ct[icell].m_s[cp].tot1 = dum1; - else + if (ct[icell].dl_s > 1e-8) { cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell); if (s_ptr) { - dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 - dum1 + min_mol); + dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 + ct[icell].m_s[cp].tot_stag - dum1 + min_mol); } } sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1; + if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol) + ct[icell].m_s[cp].tot1 = dum1 - ct[icell].m_s[cp].tot_stag - min_mol; } // check for negative moles, in the other cell... ct[icell].m_s[cp].tot2 = ct[icell].m_s[cp].tot1; - if (-ct[icell].m_s[cp].tot2 > dum2 && - (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))) + dum = 0; + if (icell == c && sptr_stag && ct[c1].m_s[cp].tot_stag) + dum = ct[c1].m_s[cp].tot_stag; + if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol && + (dV_dcell || (icell >= 0 && icell <= ilast)/* || (icell == ilast && bcon_last == 2)*/)) { - if (!ct[icell + 1].dl_s) - ct[icell].m_s[cp].tot2 = -dum2; - else + dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name); + if (ct[icell + 1].dl_s > 1e-8) { cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1); if (s_ptr) { dum2 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(ct[icell].m_s[cp].tot2 + dum2 - min_mol)); } - sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2; - if (-ct[icell].m_s[cp].tot2 > dum2) - ct[icell].m_s[cp].tot2 = -dum2; } + sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2; + if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol) + ct[icell].m_s[cp].tot2 = -dum2 + min_mol + dum; } if (fabs(ct[icell].m_s[cp].tot2) < fabs(ct[icell].m_s[cp].tot1)) ct[icell].m_s[cp].tot1 = ct[icell].m_s[cp].tot2; - - if (!cp) - { - ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0; - if (sptr_stag) - ct[i1].J_ij_sum = 0.0; - } - + if (dV_dcell || (icell > 0 && icell <= ilast)) { - dum1 -= ct[icell].m_s[cp].tot1; - if (sptr_stag) - dum1 -= ct[icell].m_s[cp].tot_stag; - sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 0e-16); - if (dum1 < -min_mol && incr > 0) + dum = ct[icell].m_s[cp].tot1; + if (stagnant) + dum += ct[icell].m_s[cp].tot_stag; + dum1 -= dum; + sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : min_mol); + if (dum1 < 0) { + dum += dum1 - min_mol; + if ((it1 = neg_moles.find(icell)) != neg_moles.end() + && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) + dum1 += it2->second; els.clear(); els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1)); neg_moles.erase(icell); neg_moles.insert(std::make_pair(icell, els)); } - else - ct[icell].J_ij_sum += dum1 * ct[icell].m_s[cp].charge; + //ct[icell].J_ij_sum -= dum * ct[icell].m_s[cp].charge; } - - if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + + if (dV_dcell || (icell >= 0 && icell < ilast)/* || (icell == ilast && bcon_last == 2)*/) { - dum2 += ct[icell].m_s[cp].tot1; - sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 0e-16); - if (dum2 < -min_mol && incr < 0) + dum = ct[icell].m_s[cp].tot1; + if (stagnant && icell == c && dV_dcell) + dum -= ct[c1].m_s[cp].tot_stag; + dum2 += dum; + sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol); + if (dum2 < 0) { + dum -= dum2 - min_mol; + if ((it1 = neg_moles.find(icell + 1)) != neg_moles.end() + && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) + dum2 += it2->second; els.clear(); els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2)); neg_moles.erase(icell + 1); neg_moles.insert(std::make_pair(icell + 1, els)); } - else - ct[icell + 1].J_ij_sum += dum2 * ct[icell].m_s[cp].charge; + //ct[icell + 1].J_ij_sum += dum * ct[icell].m_s[cp].charge; } - + if (sptr_stag) { - if ((it1 = neg_moles.find(i1)) != neg_moles.end() - && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) - { - ct[icell].m_s[cp].tot_stag += it2->second; - // the negative moles are canceled... - neg_moles.erase(it1); - els.erase(it2); - neg_moles.insert(std::make_pair(i1, els)); - } - dummy = ct[icell].m_s[cp].tot_stag; + dum = ct[icell].m_s[cp].tot_stag; if (icell == c) + dum += ct[c1].m_s[cp].tot_stag; + if (dum_stag + dum < 0) { - if (dV_dcell) - { - if ((dum = sptr2->Get_totals()[ct[icell + 1].m_s[cp].name] - ct[icell + 1].m_s[cp].tot_stag) > 0) - { - sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum; - dummy += ct[icell + 1].m_s[cp].tot_stag; - } - else - { - dummy += sptr2->Get_totals()[ct[icell].m_s[cp].name]; - sptr2->Get_totals()[ct[icell].m_s[cp].name] = 0.0; - if (dum < -min_mol) - { - els.clear(); - els.insert(std::make_pair(ct[icell].m_s[cp].name, dum)); - neg_moles.erase(icell); - neg_moles.insert(std::make_pair(icell, els)); - } - } - } - else - dummy += ct[icell + 1].m_s[cp].tot_stag; - } - if (-dummy > dum_stag) - { + dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name); if (ct[i1].dl_s) { cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, i1); if (s_ptr) { - dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dummy + dum_stag - min_mol)); + dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dum + dum_stag - min_mol)); } sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag; } } - dum_stag += dummy; - sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : 0e-16); - if (dum_stag < -min_mol) + dum_stag += dum; + sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : min_mol); + if (dum_stag < 0) { + dum -= dum_stag - min_mol; + if ((it1 = neg_moles.find(i1)) != neg_moles.end() + && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) + dum_stag += it2->second; els.clear(); els.insert(std::make_pair(ct[icell].m_s[cp].name, dum_stag)); neg_moles.erase(i1); neg_moles.insert(std::make_pair(i1, els)); } - else - ct[i1].J_ij_sum += dum_stag * ct[icell].m_s[cp].charge; + //ct[i1].J_ij_sum += dum * ct[icell].m_s[cp].charge; } - //if (cp == count_m_s - 1) + + // reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted... + if (icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1) + { + dummy = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_totals()[ct[0].m_s[cp].name] / ct[0].kgw * (1 - ct[0].dl_s); + if (dummy > 1e-6) + { + sptr1 = Utilities::Rxn_find(Rxn_solution_map, 1); + sptr2 = Utilities::Rxn_find(Rxn_solution_map, 2); + dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[1].kgw * (1 - ct[1].dl_s) - dummy; + dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[2].kgw * (1 - ct[2].dl_s) - dummy; + if (dum1 / dum2 < 0 || dum1 / dum2 > 1) + { + dum = cell_data[1].mid_cell_x / cell_data[2].mid_cell_x; + //ct[1].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[1].m_s[cp].charge; + dum1 = (dummy + dum * dum2) * ct[1].kgw / (1 - ct[1].dl_s); + sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1; + //ct[1].J_ij_sum += dum1 * ct[1].m_s[cp].charge; + } + } + dummy = Utilities::Rxn_find(Rxn_solution_map, c1)->Get_totals()[ct[0].m_s[cp].name] / ct[c1].kgw * (1 - ct[c1].dl_s); + if (dummy > 1e-6) + { + sptr1 = Utilities::Rxn_find(Rxn_solution_map, c); + sptr2 = Utilities::Rxn_find(Rxn_solution_map, c_1); + dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[c].kgw * (1 - ct[c].dl_s) - dummy; + dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[c_1].kgw * (1 - ct[c_1].dl_s) - dummy; + if (dum1 / dum2 < 0 || dum1 / dum2 > 1) + { + dum = (cell_data[c].mid_cell_x - cell_data[c_1].mid_cell_x) / + (cell_data[c1].mid_cell_x - cell_data[c_1].mid_cell_x); + //ct[c].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[c].m_s[cp].charge; + dum1 = (dummy + (1 - dum) * dum2) * ct[c].kgw / (1 - ct[c].dl_s); + sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1; + //ct[c].J_ij_sum += dum1 * ct[c].m_s[cp].charge; + } + } + } +//if (cp == count_m_s - 1) //{ // sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum); // sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum); @@ -3664,6 +3657,46 @@ 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) @@ -3848,7 +3881,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) if (dl_aq1 > 0) ct[icell].dl_s = dl_aq1 / t_aq1; if (dl_aq2 > 0) - ct[icell].dl_s = ct[jcell].dl_s = dl_aq2 / t_aq2; + { + ct[jcell].dl_s = dl_aq2 / t_aq2; + if (!ct[icell].dl_s) + ct[icell].dl_s = 1e-8; // used in implicit appt + } if (il_calcs) { @@ -4142,7 +4179,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() / mass_water_bulk_x; + dum1 = it_sc->Get_mass_water() / t_aq2; dum2 = it_sc->Get_z_gMCD_map()[1] / dum1; g_j += pow(dum2, ct[icell].v_m[k].z) * dum1; } @@ -4152,32 +4189,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } } - 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 */ + 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; else { - 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 */ + 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 (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; + calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); k++; } @@ -4249,7 +4272,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() / mass_water_bulk_x; + dum1 = it_sc->Get_mass_water() / t_aq1; dum2 = it_sc->Get_z_gMCD_map()[1] / dum1; g_i += pow(dum2, ct[icell].v_m[k].z) * dum1; } @@ -4266,31 +4289,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) g_j *= sol_D[jcell].spec[j].erm_ddl; } } - 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 */ + 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; else { - 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 */ + 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 (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; + calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant); k++; } @@ -4373,28 +4383,9 @@ 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 * (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; + 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); //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)