diff --git a/phreeqcpp/PBasic.cpp b/phreeqcpp/PBasic.cpp index f87f2a6c..09016cd9 100644 --- a/phreeqcpp/PBasic.cpp +++ b/phreeqcpp/PBasic.cpp @@ -1418,6 +1418,12 @@ listtokens(FILE * f, tokenrec * l_buf) case tokm0: output_msg("M0"); break; + case tokmcd_jtot: + output_msg("MCD_JTOT"); + break; + case tokmcd_jconc: + output_msg("MCD_JCONC"); + break; case tokmisc1: output_msg("MISC1"); break; @@ -2598,6 +2604,29 @@ factor(struct LOC_exec * LINK) } break; + case tokmcd_jtot: + { + double f = 0.0; + const char* str = stringfactor(STR1, LINK); + if (PhreeqcPtr->state == TRANSPORT && PhreeqcPtr->multi_Dflag) + { + f = PhreeqcPtr->flux_mcd(str, 1); + } + n.UU.val = (parse_all) ? 1 : f; + } + break; + case tokmcd_jconc: + { + double f = 0.0; + const char* str = stringfactor(STR1, LINK); + if (PhreeqcPtr->state == TRANSPORT && PhreeqcPtr->multi_Dflag) + { + f = PhreeqcPtr->flux_mcd(str, 2); + } + n.UU.val = (parse_all) ? 1 : f; + } + break; + case tokgamma: { const char* str = stringfactor(STR1, LINK); @@ -7436,6 +7465,8 @@ const std::map::value_type temp_tokens[] std::map::value_type("ltrim", PBasic::tokltrim), std::map::value_type("m0", PBasic::tokm0), std::map::value_type("m", PBasic::tokm), + std::map::value_type("mcd_jtot", PBasic::tokmcd_jtot), + std::map::value_type("mcd_jconc", PBasic::tokmcd_jconc), std::map::value_type("misc1", PBasic::tokmisc1), std::map::value_type("misc2", PBasic::tokmisc2), std::map::value_type("mol", PBasic::tokmol), diff --git a/phreeqcpp/PBasic.h b/phreeqcpp/PBasic.h index b3708b37..749b7077 100644 --- a/phreeqcpp/PBasic.h +++ b/phreeqcpp/PBasic.h @@ -278,6 +278,8 @@ public: tokltrim, tokm, tokm0, + tokmcd_jtot, + tokmcd_jconc, tokmisc1, tokmisc2, tokmol, diff --git a/phreeqcpp/PHRQ_io_output.cpp b/phreeqcpp/PHRQ_io_output.cpp index 46f66896..5eaeb5c1 100644 --- a/phreeqcpp/PHRQ_io_output.cpp +++ b/phreeqcpp/PHRQ_io_output.cpp @@ -201,7 +201,9 @@ void Phreeqc:: screen_msg(const char *err_str) /* ---------------------------------------------------------------------- */ { +#ifndef TESTING if (phrq_io) phrq_io->screen_msg(err_str); +#endif } // ---------------------------------------------------------------------- */ // dump file methods diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 4023237b..7e4e0c7d 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -111,6 +111,7 @@ public: LDBLE phase_vm(const char* phase_name); LDBLE diff_c(const char* species_name); LDBLE setdiff_c(const char* species_name, double d); + LDBLE flux_mcd(const char* species_name, int option); LDBLE sa_declercq(double type, double sa, double d, double m, double m0, double gfw); LDBLE calc_SC(void); /* VP: Density Start */ diff --git a/phreeqcpp/class_main.cpp b/phreeqcpp/class_main.cpp index ee43e385..2730fc9d 100644 --- a/phreeqcpp/class_main.cpp +++ b/phreeqcpp/class_main.cpp @@ -283,6 +283,9 @@ int Phreeqc:: write_banner(void) /* ---------------------------------------------------------------------- */ { +#ifdef TESTING + return OK; +#endif char buffer[80]; int len, indent; screen_msg( diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index b217d965..619c75e6 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -1535,6 +1535,18 @@ public: const char* name; LDBLE tot1, tot2, tot_stag, charge; }; +class J_ij_save +{ +public: + ~J_ij_save() {}; + J_ij_save() + { + // species change in cells i and j + flux_t = 0; + flux_c = 0; + } + double flux_t, flux_c; +}; class M_S { public: diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 00fedd8a..2f42f061 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -44,6 +44,7 @@ struct CT /* summed parts of V_M and mcd transfer in a timestep for all cells, f class M_S *m_s; int v_m_size, J_ij_size, m_s_size; } *ct = NULL; +std::map > cell_J_ij; struct MOLES_ADDED /* total moles added to balance negative conc's */ { char *name; @@ -501,7 +502,7 @@ transport(void) else sprintf(token, "\nCalculating transport: %d (mobile) cells, %d shifts, %d mixruns...\n\n", count_cells, count_shifts - transport_start + 1, nmix); - screen_msg(token); + warning_msg(token); max_iter = 0; for (transport_step = transport_start; transport_step <= count_shifts; transport_step++) @@ -960,8 +961,7 @@ transport(void) { sprintf(token, "\nFor balancing negative concentrations in MCD, added in total to the system:"); - if (phrq_io) - phrq_io->warning_msg(token); + warning_msg(token); for (i = 0; i < count_moles_added; i++) { if (!moles_added[i].moles) @@ -969,8 +969,7 @@ transport(void) sprintf(token, "\t %.4e moles %s.", (double)moles_added[i].moles, moles_added[i].name); - if (phrq_io) - phrq_io->warning_msg(token); + warning_msg(token); } } } @@ -1048,7 +1047,7 @@ transport_cleanup(void) mixf[i] = (LDBLE *)free_check_null(mixf[i]); if (l_stag) mixf_stag[i] = (LDBLE *)free_check_null(mixf_stag[i]); - if (!dV_dcell) + if (!dV_dcell && !fix_current) { cell_data[i].potV = 0.0; use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i)); @@ -1078,7 +1077,7 @@ print_punch(int i, boolean active) if (!active) run_reactions(i, 0, NOMIX, 0); cell_no = i; - if (dV_dcell) + if (dV_dcell || fix_current) { use.Set_n_solution_user(i); use.Get_solution_ptr()->Set_potV(cell_data[i].potV); @@ -2238,6 +2237,17 @@ diffuse_implicit(LDBLE DDt, int stagnant) int c = count_cells, c1 = c + 1, c2 = c + 2, c_1 = c - 1, cc = c + stagnant * c, cc1 = cc + 1; comp = sol_D[1].count_spec - sol_D[1].count_exch_spec; + cell_J_ij.clear(); + for (i = 0; i <= count_cells; i++) + { + std::map J_map; + for (cp = 0; cp < comp; cp++) + { + class J_ij_save J_save; + J_map[sol_D[1].spec[cp].name] = J_save; + } + cell_J_ij[i] = J_map; + } if (heat_nmix) comp += 1; @@ -2733,7 +2743,16 @@ diffuse_implicit(LDBLE DDt, int stagnant) ct[i].J_ij[cp].charge = ct[i].v_m[cp].z; grad = (Ct2[i1] - Ct2[i])/* * ct[i].v_m[cp].D*/; // .D has the d{lg(gamma)} / d{lg(molal)} correction if (!i0) + { ct[i].J_ij[cp].tot1 = -mixf[i][cp] * grad; + + std::map >::iterator + cell_iter = cell_J_ij.find(i); + std::map< std::string, class J_ij_save >::iterator + s_iter = cell_iter->second.find(sol_D[1].spec[cp].name); + s_iter->second.flux_c = ct[i].J_ij[cp].tot1; + s_iter->second.flux_t = ct[i].J_ij[cp].tot1; + } else ct[i].J_ij[cp].tot_stag = -mixf_stag[i][cp] * grad; if (ct[i].v_m[cp].z) @@ -2796,7 +2815,7 @@ diffuse_implicit(LDBLE DDt, int stagnant) sum_Rd += (current_cells[0].dif - current_cells[i].dif) * current_cells[i].R; } } - if (dV_dcell) + if (dV_dcell || fix_current) { if (fix_current) { @@ -2848,6 +2867,12 @@ diffuse_implicit(LDBLE DDt, int stagnant) { dVc = (cell_data[i + 1].potV - cell_data[i].potV) * F_Re3 / l_tk_x2[i]; ct[i].J_ij[cp].tot1 -= mixf[i][cp] * dVc; + + std::map >::iterator + cell_iter = cell_J_ij.find(i); + std::map< std::string, class J_ij_save >::iterator + s_iter = cell_iter->second.find(sol_D[1].spec[cp].name); + s_iter->second.flux_t = ct[i].J_ij[cp].tot1; } if (stagnant && ct[i].Dz2c_stag) { @@ -3281,7 +3306,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) dVtemp = dV_dcell; dV_dcell = 0; } - + cell_J_ij.clear(); icell = jcell = -1; first_c = last_c = -1; il_calcs = -1; @@ -4485,21 +4510,35 @@ dV_dcell2: if (ct[icell].v_m[i].z) Sum_zM += ct[icell].v_m[i].b_ij * ct[icell].v_m[i].z * ct[icell].v_m[i].grad; } + std::map J_map; for (i = 0; i < ct[icell].J_ij_count_spec; i++) { + class J_ij_save J_save; ct[icell].J_ij[i].tot1 = -ct[icell].v_m[i].grad; + J_save.flux_c = ct[icell].J_ij[i].tot1; ct[icell].J_ij[i].charge = ct[icell].v_m[i].z; if (!dV_dcell && ct[icell].v_m[i].z && ct[icell].Dz2c > 0) { ct[icell].J_ij[i].tot1 += Sum_zM * ct[icell].v_m[i].zc / ct[icell].Dz2c; } + J_save.flux_t = ct[icell].J_ij[i].tot1; if (stagnant) + { ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * 2 * mixf; + J_save.flux_c *= ct[icell].v_m[i].b_ij * 2 * mixf; + J_save.flux_t *= ct[icell].v_m[i].b_ij * 2 * mixf; + } else + { ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * DDt; + J_save.flux_c *= ct[icell].v_m[i].b_ij * DDt; + J_save.flux_t *= ct[icell].v_m[i].b_ij * DDt; + } + J_map[ct[icell].J_ij[i].name] = J_save; ct[icell].J_ij[i].tot2 = ct[icell].J_ij[i].tot1; //ct[icell].J_ij_sum += ct[icell].v_m[i].z * ct[icell].J_ij[i].tot1; } + cell_J_ij[icell] = J_map; // assure that icell and jcell have dl water when checking negative conc's in MCD ct[icell].dl_s = dl_aq1; ct[jcell].dl_s = dl_aq2; @@ -6104,3 +6143,38 @@ calc_vm_Cl(void) } return V_Cl; } +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +flux_mcd(const char* species_name, int option) +/* ---------------------------------------------------------------------- */ +{ + class species* s_ptr; + double f = 0.0, dum = 0.0; + if (state == TRANSPORT && multi_Dflag) + { + s_ptr = s_search(species_name); + if (s_ptr != NULL && s_ptr->in != FALSE && s_ptr->type < EMINUS) + { + int n = cell_no; + std::map >::iterator + j_map_iter = cell_J_ij.find(n); + if (j_map_iter != cell_J_ij.end()) + { + std::map::iterator + j_save_iter = j_map_iter->second.find(species_name); + if (j_save_iter != j_map_iter->second.end()) + { + if (option == 1) + { + f = j_save_iter->second.flux_t; + } + else if (option == 2) + { + f = j_save_iter->second.flux_c; + } + } + } + } + } + return (f); +} \ No newline at end of file