diff --git a/phreeqcpp/ExchComp.cxx b/phreeqcpp/ExchComp.cxx index 5af520b7..dba46b3d 100644 --- a/phreeqcpp/ExchComp.cxx +++ b/phreeqcpp/ExchComp.cxx @@ -360,7 +360,8 @@ cxxExchComp::add(const cxxExchComp & addee, LDBLE extensive) this->la = f1 * this->la + f2 * addee.la; this->charge_balance += addee.charge_balance * extensive; - if (this->phase_name != addee.phase_name) + //if (this->phase_name != addee.phase_name) + if (Phreeqc::strcmp_nocase(this->phase_name.c_str(), addee.phase_name.c_str()) != 0) { std::ostringstream oss; oss << @@ -374,7 +375,8 @@ cxxExchComp::add(const cxxExchComp & addee, LDBLE extensive) this->phase_proportion = this->phase_proportion * f1 + addee.phase_proportion * f2; } - if (this->rate_name != addee.rate_name) + //if (this->rate_name != addee.rate_name) + if (Phreeqc::strcmp_nocase(this->rate_name.c_str(), addee.rate_name.c_str()) != 0) { std::ostringstream oss; oss << diff --git a/phreeqcpp/Exchange.cxx b/phreeqcpp/Exchange.cxx index 117123dd..c1192df9 100644 --- a/phreeqcpp/Exchange.cxx +++ b/phreeqcpp/Exchange.cxx @@ -467,7 +467,7 @@ cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std: this->exchange_comps.clear(); for (int n = 0; n < count; n++) { - cxxExchComp ec; + cxxExchComp ec(this->io); ec.Deserialize(dictionary, ints, doubles, ii, dd); this->exchange_comps.push_back(ec); } diff --git a/phreeqcpp/PPassemblage.cxx b/phreeqcpp/PPassemblage.cxx index 9af480db..3536dc51 100644 --- a/phreeqcpp/PPassemblage.cxx +++ b/phreeqcpp/PPassemblage.cxx @@ -373,7 +373,7 @@ cxxPPassemblage::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->pp_assemblage_comps.clear(); for (int n = 0; n < count; n++) { - cxxPPassemblageComp ppc; + cxxPPassemblageComp ppc(this->io); ppc.Deserialize(dictionary, ints, doubles, ii, dd); std::string str(ppc.Get_name()); this->pp_assemblage_comps[str] = ppc; diff --git a/phreeqcpp/Phreeqc.cpp b/phreeqcpp/Phreeqc.cpp index a2b9bdaf..e300d820 100644 --- a/phreeqcpp/Phreeqc.cpp +++ b/phreeqcpp/Phreeqc.cpp @@ -713,6 +713,7 @@ void Phreeqc::init(void) interlayer_Dflag = FALSE; implicit = FALSE; max_mixf = 1.0; + min_dif_LM = -30.0; default_Dw = 0; correct_Dw = 0; multi_Dpor = 0; @@ -1691,6 +1692,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) interlayer_Dflag = pSrc->interlayer_Dflag; implicit = pSrc->implicit; max_mixf = pSrc->max_mixf; + min_dif_LM = pSrc->min_dif_LM; default_Dw = pSrc->default_Dw; correct_Dw = pSrc->correct_Dw; multi_Dpor = pSrc->multi_Dpor; diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index 96df04da..21ef40b5 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -1066,9 +1066,11 @@ 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 diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant); - int fill_spec(int cell_no); - void define_ct_structures(void); + 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); + LDBLE moles_from_donnan_layer(cxxSurface *sptr, const char *name, LDBLE moles_needed); + LDBLE add_MCD_moles(LDBLE moles, LDBLE min_mol, int i, cxxSolution *sptr, const char *name); int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i, int stagnant); static int sort_species_name(const void *ptr1, const void *ptr2); int disp_surf(LDBLE stagkin_time); @@ -1429,6 +1431,7 @@ protected: int interlayer_Dflag; /* multicomponent diffusion and diffusion through interlayer porosity */ int implicit; /* implicit calculation of diffusion */ LDBLE max_mixf; /* the maximum value of the implicit mixfactor = De * Dt / (Dx^2) */ + LDBLE min_dif_LM; /* the minimal log10(molality) for including a species in multicomponent diffusion */ LDBLE default_Dw; /* default species diffusion coefficient in water at 25oC, m2/s */ int correct_Dw; /* if true, Dw is adapted in calc_SC */ LDBLE multi_Dpor; /* uniform porosity of free porewater in solid medium */ diff --git a/phreeqcpp/SS.cxx b/phreeqcpp/SS.cxx index 2ed2edc6..a40e577f 100644 --- a/phreeqcpp/SS.cxx +++ b/phreeqcpp/SS.cxx @@ -580,7 +580,7 @@ cxxSS::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->ss_comps.clear(); for (int i = 0; i < count; i++) { - cxxSScomp ssc; + cxxSScomp ssc(this->io); ssc.Deserialize(dictionary, ints, doubles, ii, dd); this->ss_comps.push_back(ssc); } diff --git a/phreeqcpp/SSassemblage.cxx b/phreeqcpp/SSassemblage.cxx index 471127bb..ef945efa 100644 --- a/phreeqcpp/SSassemblage.cxx +++ b/phreeqcpp/SSassemblage.cxx @@ -319,7 +319,7 @@ cxxSSassemblage::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->SSs.clear(); for (int n = 0; n < count; n++) { - cxxSS ssc; + cxxSS ssc(this->io); ssc.Deserialize(dictionary, ints, doubles, ii, dd); std::string str(ssc.Get_name()); this->SSs[str] = ssc; diff --git a/phreeqcpp/Serializer.cxx b/phreeqcpp/Serializer.cxx index 4b67baeb..7b05d817 100644 --- a/phreeqcpp/Serializer.cxx +++ b/phreeqcpp/Serializer.cxx @@ -133,7 +133,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto break; case PT_GASPHASE: { - cxxGasPhase entity; + cxxGasPhase entity(phreeqc_ref.Get_phrq_io()); entity.Deserialize(dictionary, ints, doubles, ii, dd); int n_user = entity.Get_n_user(); phreeqc_ref.Get_Rxn_gas_phase_map()[n_user] = entity; @@ -141,7 +141,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto break; case PT_KINETICS: { - cxxKinetics entity; + cxxKinetics entity(phreeqc_ref.Get_phrq_io()); entity.Deserialize(dictionary, ints, doubles, ii, dd); int n_user = entity.Get_n_user(); phreeqc_ref.Get_Rxn_kinetics_map()[n_user] = entity; @@ -158,7 +158,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto break; case PT_SSASSEMBLAGE: { - cxxSSassemblage entity; + cxxSSassemblage entity(phreeqc_ref.Get_phrq_io()); entity.Deserialize(dictionary, ints, doubles, ii, dd); int n_user = entity.Get_n_user(); phreeqc_ref.Get_Rxn_ss_assemblage_map()[n_user] = entity; diff --git a/phreeqcpp/StorageBin.cxx b/phreeqcpp/StorageBin.cxx index 24e698d0..cf4d5c56 100644 --- a/phreeqcpp/StorageBin.cxx +++ b/phreeqcpp/StorageBin.cxx @@ -1061,7 +1061,7 @@ cxxStorageBin::read_raw(CParser & parser) break; case Keywords::KEY_SOLID_SOLUTIONS_RAW: { - cxxSSassemblage entity; + cxxSSassemblage entity(this->Get_io()); entity.read_raw(parser); SSassemblages[entity.Get_n_user()] = entity; } @@ -1195,7 +1195,7 @@ cxxStorageBin::read_raw_keyword(CParser & parser) case Keywords::KEY_SOLID_SOLUTIONS_RAW: { - cxxSSassemblage entity; + cxxSSassemblage entity(this->Get_io()); entity.read_raw(parser); SSassemblages[entity.Get_n_user()] = entity; entity_number = entity.Get_n_user(); diff --git a/phreeqcpp/Surface.cxx b/phreeqcpp/Surface.cxx index 99a986de..fdbc122a 100644 --- a/phreeqcpp/Surface.cxx +++ b/phreeqcpp/Surface.cxx @@ -812,7 +812,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->surface_comps.clear(); for (int n = 0; n < count; n++) { - cxxSurfaceComp sc; + cxxSurfaceComp sc(this->io); sc.Deserialize(dictionary, ints, doubles, ii, dd); this->surface_comps.push_back(sc); } @@ -822,7 +822,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->surface_charges.clear(); for (int n = 0; n < count; n++) { - cxxSurfaceCharge sc; + cxxSurfaceCharge sc(this->io); sc.Deserialize(dictionary, ints, doubles, ii, dd); this->surface_charges.push_back(sc); } diff --git a/phreeqcpp/SurfaceComp.cxx b/phreeqcpp/SurfaceComp.cxx index 5a23eb16..8ef1669c 100644 --- a/phreeqcpp/SurfaceComp.cxx +++ b/phreeqcpp/SurfaceComp.cxx @@ -395,7 +395,8 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive) this->totals.add_extensive(addee.totals, extensive); this->la = f1 * this->la + f2 * addee.la; this->charge_balance += addee.charge_balance * extensive; - if (this->phase_name != addee.phase_name) + //if (this->phase_name != addee.phase_name) + if (Phreeqc::strcmp_nocase(this->phase_name.c_str(),addee.phase_name.c_str()) != 0) { std::ostringstream oss; oss << @@ -410,11 +411,12 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive) this->phase_proportion * f1 + addee.phase_proportion * f2; } - if (this->rate_name != addee.rate_name) + //if (this->rate_name != addee.rate_name) + if (Phreeqc::strcmp_nocase(this->rate_name.c_str(), addee.rate_name.c_str()) != 0) { std::ostringstream oss; oss << - "Cannot mix two exchange components with same formula and different related kinetics, " + "Cannot mix two surface components with same formula and different related kinetics, " << this->formula; error_msg(oss.str().c_str(), CONTINUE); return; @@ -430,7 +432,7 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive) { std::ostringstream oss; oss << - "Cannot mix exchange components related to phase with exchange components related to kinetics, " + "Cannot mix surface components related to phase with exchange components related to kinetics, " << this->formula; error_msg(oss.str().c_str(), CONTINUE); return; diff --git a/phreeqcpp/class_main.cpp b/phreeqcpp/class_main.cpp index 9ad59d8f..2ee57fc6 100644 --- a/phreeqcpp/class_main.cpp +++ b/phreeqcpp/class_main.cpp @@ -277,8 +277,7 @@ write_banner(void) /* version */ #ifdef NPP - //len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1 AmpŠre"); - len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1"); + len = sprintf(buffer, "* PHREEQC-%s *", "3.5.2"); #else len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@"); #endif @@ -302,7 +301,7 @@ write_banner(void) /* date */ #ifdef NPP - len = sprintf(buffer, "%s", "May 29, 2019"); + len = sprintf(buffer, "%s", "August 1, 2019"); #else len = sprintf(buffer, "%s", "@VER_DATE@"); #endif @@ -491,13 +490,10 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, user_database = string_duplicate(token); screen_msg(sformatf("Database file: %s\n\n", token)); strcpy(db_file, token); -#ifdef NPP - //output_msg(sformatf("Using PHREEQC: version 3.5.1 Ampère, compiled May 29, 2019\n")); -#endif output_msg(sformatf(" Input file: %s\n", in_file)); output_msg(sformatf(" Output file: %s\n", out_file)); #ifdef NPP - output_msg(sformatf("Using PHREEQC: version 3.5.1, compiled May 29, 2019\n")); + output_msg(sformatf("Using PHREEQC: version 3.5.2, compiled August 1, 2019\n")); #endif output_msg(sformatf("Database file: %s\n\n", token)); #ifdef NPP diff --git a/phreeqcpp/common/PHRQ_io.cpp b/phreeqcpp/common/PHRQ_io.cpp index c57c4b7f..876a92f0 100644 --- a/phreeqcpp/common/PHRQ_io.cpp +++ b/phreeqcpp/common/PHRQ_io.cpp @@ -792,6 +792,7 @@ get_line(void) continue; #else output_msg(errstr.str().c_str()); + output_flush(); error_msg(errstr.str().c_str(), OT_STOP); #endif } diff --git a/phreeqcpp/cvode.cpp b/phreeqcpp/cvode.cpp index 171ce1c1..58c5bccd 100644 --- a/phreeqcpp/cvode.cpp +++ b/phreeqcpp/cvode.cpp @@ -72,6 +72,7 @@ #include "sundialsmath.h" #include "Phreeqc.h" + #if !defined(WIN32_MEMORY_DEBUG) #define malloc MACHENV_MALLOC PHRQ_malloc #endif @@ -1974,7 +1975,6 @@ CVStep(CVodeMem cv_mem) } #endif CVSet(cv_mem); - nflag = CVnls(cv_mem, nflag); if (CVMEM cvode_error == TRUE || predict_fail) { @@ -2785,7 +2785,6 @@ CVnlsNewton(CVodeMem cv_mem, int nflag) loop { - f(N, tn, zn[0], ftemp, f_data); nfe++; @@ -3408,7 +3407,7 @@ CVHandleFailure(CVodeMem cv_mem, int kflag) size is reset accordingly. *****************************************************************/ - +#ifdef ORIGINAL void CVBDFStab(CVodeMem cv_mem) { @@ -3469,6 +3468,82 @@ CVBDFStab(CVodeMem cv_mem) nscon = 0; } } +#endif +void +CVBDFStab(CVodeMem cv_mem) +{ + // appt try... + if (q >= 3 && qprime >= q) + { + if (tq[5] < saved_tq5) + qprime = 1; + //else + //nscon = 0; + } + + //int i, k, ldflag, factorial; + //realtype sq, sqm1, sqm2; + + ///* If order is 3 or greater, then save scaled derivative data, + // push old data down in i, then add current values to top. */ + + //if (q >= 3) + //{ + // for (k = 1; k <= 3; k++) + // { + // for (i = 5; i >= 2; i--) + // ssdat[i][k] = ssdat[i - 1][k]; + // } + // factorial = 1; + // for (i = 1; i <= q - 1; i++) + // factorial *= i; + // sq = factorial * q * (q + 1) * acnrm / MAX(tq[5], TINY); + // sqm1 = factorial * q * N_VWrmsNorm(zn[q], ewt); + // sqm2 = factorial * N_VWrmsNorm(zn[q - 1], ewt); + // ssdat[1][1] = sqm2 * sqm2; + // ssdat[1][2] = sqm1 * sqm1; + // ssdat[1][3] = sq * sq; + //} + + //if (qprime >= q) + //{ + + // /* If order is 3 or greater, and enough ssdat has been saved, + // nscon >= q+5, then call stability limit detection routine. */ + + // if ((q >= 3) && (nscon >= q + 5)) + // { + // ldflag = CVsldet(cv_mem); + // //cv_mem->cv_machenv->phreeqc_ptr->set_forward_output_to_log(1); // appt + // qprime = 1; // appt try + // //CVMEM warning_msg(CVMEM sformatf( + // // "CVBDFStab: ldflag = %d, order(q) = %d, qprime = %d, nst = %d, h = %8.2e, time = %8.2e\n", + // // ldflag, q, qprime, nst, h, CVMEM cvode_last_good_time)); + + // if (ldflag > 3) + // { + // /* A stability limit violation is indicated by + // a return flag of 4, 5, or 6. + // Reduce new order. */ + // qprime = q - 1; + // eta = etaqm1; + // eta = MIN(eta, etamax); + // eta = eta / MAX(ONE, ABS(h) * hmax_inv * eta); + // hprime = h * eta; + // iopt[NOR] = iopt[NOR] + 1; + // //CVMEM warning_msg(CVMEM sformatf( + // // " Order reduced to %d by CVBDFStab at nst = %d,\n h = %e hnew = %e\n", + // // qprime,nst,h,h*eta)); + // } + // } + //} + //else + //{ + // /* Otherwise, let order increase happen, and + // reset stability limit counter, nscon. */ + // nscon = 0; + //} +} /********************* CVsldet ************************************ This routine detects stability limitation using stored scaled diff --git a/phreeqcpp/cxxKinetics.cxx b/phreeqcpp/cxxKinetics.cxx index 17ad6626..be640cae 100644 --- a/phreeqcpp/cxxKinetics.cxx +++ b/phreeqcpp/cxxKinetics.cxx @@ -605,7 +605,7 @@ cxxKinetics::Deserialize(Dictionary & dictionary, std::vector < int >&ints, this->kinetics_comps.clear(); for (int i = 0; i < n; i++) { - cxxKineticsComp kc; + cxxKineticsComp kc(this->io); kc.Deserialize(dictionary, ints, doubles, ii, dd); this->kinetics_comps.push_back(kc); } diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h index 7b759984..29be4855 100644 --- a/phreeqcpp/global_structures.h +++ b/phreeqcpp/global_structures.h @@ -547,6 +547,7 @@ struct cell_data LDBLE potV; /* potential (V) */ int punch; int print; + int same_model; }; /*---------------------------------------------------------------------- @@ -1052,12 +1053,12 @@ struct sol_D struct J_ij { const char *name; - LDBLE tot1, tot2, charge; /* species change in cells i and j */ + LDBLE tot1, tot2, tot_stag, charge; /* species change in cells i and j */ }; struct M_S { const char *name; - LDBLE tot1, tot2, charge; /* master species transport in cells i and j */ + LDBLE tot1, tot2, tot_stag, charge; /* master species transport in cells i and j */ }; // Pitzer definitions typedef enum diff --git a/phreeqcpp/integrate.cpp b/phreeqcpp/integrate.cpp index 8a7cb408..3636b1ed 100644 --- a/phreeqcpp/integrate.cpp +++ b/phreeqcpp/integrate.cpp @@ -1069,9 +1069,17 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq) l_iter++; if (l_iter > 50) { + pr.all = TRUE; + pr.exchange = TRUE; + pr.headings = TRUE; + pr.pp_assemblage = TRUE; + pr.species = TRUE; + pr.surface = TRUE; + pr.totals = TRUE; + print_all(); error_string = sformatf( - "\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n", - (double) surf_chrg_eq, (double) charge_ptr->Get_mass_water()); + "\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n", + (double)surf_chrg_eq, (double)charge_ptr->Get_mass_water()); error_msg(error_string, STOP); } } diff --git a/phreeqcpp/kinetics.cpp b/phreeqcpp/kinetics.cpp index c13af7d6..6396a192 100644 --- a/phreeqcpp/kinetics.cpp +++ b/phreeqcpp/kinetics.cpp @@ -500,8 +500,8 @@ rk_kinetics(int i, LDBLE kin_time, int use_mix, int nsaver, if (step_bad > kinetics_ptr->Get_bad_step_max()) { error_string = sformatf( - "Bad RK steps > %d. Please decrease (time)step or increase -bad_step_max.", - kinetics_ptr->Get_bad_step_max()); + "Bad RK steps > %d in cell %d. Please decrease (time)step or increase -bad_step_max.", + kinetics_ptr->Get_bad_step_max(), cell_no); error_msg(error_string, STOP); } @@ -1511,11 +1511,11 @@ restart: } if (j == 14) { - cxxStorageBin error_bin; + cxxStorageBin error_bin(this->Get_phrq_io()); Use2cxxStorageBin(error_bin); std::ostringstream error_input; error_bin.dump_raw(error_input, 0); - cxxStorageBin reread; + cxxStorageBin reread(this->Get_phrq_io()); std::istringstream is(error_input.str()); CParser cp(is); cp.set_echo_stream(CParser::EO_NONE); @@ -1574,7 +1574,7 @@ restart: * write to error.inp what failed to converge. */ std::ofstream error_input("error.inp"); - cxxStorageBin error_bin; + cxxStorageBin error_bin(this->Get_phrq_io()); Use2cxxStorageBin(error_bin); error_bin.dump_raw(error_input, 0); error_input.close(); diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index f0512ccf..8b72bcb9 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -1396,7 +1396,7 @@ xss_assemblage_save(int n_user) * Save ss_assemblage composition into structure ss_assemblage with user * number n_user. */ - cxxSSassemblage temp_ss_assemblage; + cxxSSassemblage temp_ss_assemblage(this->phrq_io); if (use.Get_ss_assemblage_ptr() == NULL) return (OK); @@ -1719,6 +1719,8 @@ xsurface_save(int n_user) if (x[i]->type == SURFACE) { cxxSurfaceComp *comp_ptr = temp_surface.Find_comp(x[i]->surface_comp); + if (comp_ptr == NULL) + continue; // appt in transport with different mobile and stagnant surfaces assert(comp_ptr); comp_ptr->Set_la(x[i]->master[0]->s->la); comp_ptr->Set_moles(0.); @@ -1748,6 +1750,8 @@ xsurface_save(int n_user) else if (x[i]->type == SURFACE_CB && (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM)) { cxxSurfaceCharge *charge_ptr = temp_surface.Find_charge(x[i]->surface_charge); + if (charge_ptr == NULL) + continue; // appt in transport with different mobile and stagnant surfaces assert(charge_ptr); charge_ptr->Set_charge_balance(x[i]->f); charge_ptr->Set_la_psi(x[i]->master[0]->s->la); diff --git a/phreeqcpp/phqalloc.cpp b/phreeqcpp/phqalloc.cpp index 8b7ce14b..08020dd8 100644 --- a/phreeqcpp/phqalloc.cpp +++ b/phreeqcpp/phqalloc.cpp @@ -141,8 +141,8 @@ PHRQ_calloc(size_t num, size_t size assert((s_pTail == NULL) || (s_pTail->pNext == NULL)); - p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num); - //p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt + //p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num); + p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt if (p == NULL) return NULL; diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp index b99716d5..45446268 100644 --- a/phreeqcpp/prep.cpp +++ b/phreeqcpp/prep.cpp @@ -6034,6 +6034,8 @@ check_same_model(void) last_model.force_prep = FALSE; return (FALSE); } + if (state == TRANSPORT && cell_data[cell_no].same_model) + return TRUE; /* * Check master species */ diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp index 0525b6c9..25334704 100644 --- a/phreeqcpp/print.cpp +++ b/phreeqcpp/print.cpp @@ -1353,8 +1353,9 @@ print_pp_assemblage(void) { if (x[j]->type != PP) continue; - //cxxPPassemblageComp * comp_ptr = pp_assemblage_ptr->Find(x[j]->pp_assemblage_comp_name); - cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr; + cxxPPassemblage * pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, use.Get_n_pp_assemblage_user()); + cxxPPassemblageComp * comp_ptr = pp_assemblage_ptr->Find(x[j]->pp_assemblage_comp_name); + //cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr; // appt, is sometimes lost?? /* * Print saturation index */ diff --git a/phreeqcpp/read.cpp b/phreeqcpp/read.cpp index 9d6f7043..6b7843b3 100644 --- a/phreeqcpp/read.cpp +++ b/phreeqcpp/read.cpp @@ -1007,7 +1007,7 @@ read_exchange(void) input_error++; break; } - cxxExchComp temp_comp; + cxxExchComp temp_comp(this->phrq_io); temp_exchange.Get_exchange_comps().push_back(temp_comp); comp_ptr = &(temp_exchange.Get_exchange_comps().back()); comp_ptr->Set_formula(token.c_str()); @@ -1268,7 +1268,7 @@ read_gas_phase(void) char *ptr; char *description; char token[MAX_LENGTH]; - cxxGasPhase temp_gas_phase; + cxxGasPhase temp_gas_phase(this->phrq_io); int return_value, opt; char *next_char; const char *opt_list[] = { @@ -2057,7 +2057,7 @@ read_kinetics(void) */ ptr = line; read_number_description(ptr, &n_user, &n_user_end, &description); - cxxKinetics temp_kinetics; + cxxKinetics temp_kinetics(this->phrq_io); temp_kinetics.Set_n_user(n_user); temp_kinetics.Set_n_user_end(n_user_end); temp_kinetics.Set_description(description); @@ -7538,7 +7538,8 @@ read_surface(void) break; } - cxxSurfaceComp temp_comp; + cxxSurfaceComp temp_comp(this->phrq_io); + temp_surface.Get_surface_comps().push_back(temp_comp); comp_ptr = &(temp_surface.Get_surface_comps().back()); comp_ptr->Set_formula(token.c_str()); @@ -7648,7 +7649,7 @@ read_surface(void) formula = (char*)free_check_null(formula); if (charge_ptr == NULL) { - cxxSurfaceCharge temp_charge; + cxxSurfaceCharge temp_charge(this->phrq_io); temp_charge.Set_name(name); if (comp_ptr->Get_phase_name().size() == 0 && comp_ptr->Get_rate_name().size() == 0) @@ -10125,7 +10126,7 @@ read_solid_solutions(void) case 0: /* component */ case 1: /* comp */ { - cxxSScomp comp; + cxxSScomp comp(this->phrq_io); /* * Read phase name of component */ diff --git a/phreeqcpp/readtr.cpp b/phreeqcpp/readtr.cpp index fee12eea..897edf78 100644 --- a/phreeqcpp/readtr.cpp +++ b/phreeqcpp/readtr.cpp @@ -37,13 +37,13 @@ read_transport(void) */ char *ptr; int i, j, l; - int count_length, count_disp, count_punch, count_print, count_por; + int count_length, count_disp, count_punch, count_print, count_por, count_same_model; int count_length_alloc, count_disp_alloc, count_por_alloc; char token[MAX_LENGTH]; char *description; int n_user, n_user_end; LDBLE *length, *disp, *pors; - int *punch_temp, *print_temp; + int *punch_temp, *print_temp, *same_model_temp; int return_value, opt, opt_save; char *next_char, *next_char_save; char file_name[MAX_LENGTH]; @@ -95,9 +95,10 @@ read_transport(void) "porosity", /* 43 */ "fix_current", /* 44 */ "current", /* 45 */ - "implicit" /* 46 */ + "implicit", /* 46 */ + "same_model" /* 47 */ }; - int count_opt_list = 47; + int count_opt_list = 48; strcpy(file_name, "phreeqc.dmp"); /* @@ -113,7 +114,7 @@ read_transport(void) } else old_cells = count_cells; - count_length = count_disp = count_punch = count_print = count_por = 0; + count_length = count_disp = count_punch = count_print = count_por = count_same_model = 0; length = (LDBLE *)PHRQ_malloc(sizeof(LDBLE)); if (length == NULL) @@ -135,6 +136,10 @@ read_transport(void) if (print_temp == NULL) malloc_error(); + same_model_temp = (int *)PHRQ_malloc(sizeof(int)); + if (same_model_temp == NULL) + malloc_error(); + count_length_alloc = count_disp_alloc = count_por_alloc = 1; transport_start = 1; /* @@ -706,8 +711,27 @@ read_transport(void) //warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion."); max_mixf = 1.0; } + min_dif_LM = -30.0; + if (copy_token(token, &next_char, &l) != EMPTY) + { + /* minimal moles for diffusion */ + if (sscanf(token, SCANFORMAT, &min_dif_LM) != 1) + { + input_error++; + error_string = sformatf( + "Expected the minimal log10(molality) for including a species in multicomponent diffusion,\n taking -30.0"); + warning_msg(error_string); + break; + } + } opt_save = OPTION_DEFAULT; break; + case 47: /* same_model */ + same_model_temp = + read_list_ints_range(&next_char, &count_same_model, FALSE, + same_model_temp); + opt_save = 47; + break; } if (return_value == EOF || return_value == KEYWORD) break; @@ -767,6 +791,7 @@ read_transport(void) cell_data[i].potV = 0; cell_data[i].punch = FALSE; cell_data[i].print = FALSE; + cell_data[i].same_model = FALSE; } all_cells = all_cells_now; } @@ -961,6 +986,29 @@ read_transport(void) else if (simul_tr == 1 || old_cells != count_cells) for (i = 0; i < all_cells; i++) cell_data[i].print = TRUE; + /* + * Fill in data for same_model + */ + if (count_same_model != 0) + { + for (i = 0; i < all_cells; i++) + cell_data[i].same_model = FALSE; + for (i = 0; i < count_same_model; i++) + { + if (same_model_temp[i] > all_cells - 1 || same_model_temp[i] < 0) + { + error_string = sformatf( + "Cell number for same_model is out of range, %d. Request ignored.", + same_model_temp[i]); + warning_msg(error_string); + } + else + cell_data[same_model_temp[i]].same_model = TRUE; + } + } + else if (simul_tr == 1 || old_cells != count_cells) + for (i = 0; i < all_cells; i++) + cell_data[i].same_model = FALSE; //#define OLD_POROSITY #if defined(OLD_POROSITY) /* @@ -1043,7 +1091,7 @@ read_transport(void) "Must enter diffusion coefficient (-diffc) when modeling thermal diffusion."); error_msg(error_string, CONTINUE); } - else if (heat_diffc > diffc) + else if (heat_diffc > diffc && !implicit) { error_string = sformatf( "Thermal diffusion is calculated assuming exchange factor was for\n\t effective (non-thermal) diffusion coefficient = %e.", @@ -1053,7 +1101,7 @@ read_transport(void) } else { - if (heat_diffc > diffc) + if (heat_diffc > diffc && !implicit) { input_error++; error_string = sformatf( @@ -1077,6 +1125,7 @@ read_transport(void) pors = (LDBLE *)free_check_null(pors); punch_temp = (int *)free_check_null(punch_temp); print_temp = (int *)free_check_null(print_temp); + same_model_temp = (int *)free_check_null(same_model_temp); if (dump_in == TRUE) { diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp index 12c2c7f1..3eda7432 100644 --- a/phreeqcpp/transport.cpp +++ b/phreeqcpp/transport.cpp @@ -20,7 +20,7 @@ std::set dif_spec_names; std::set dif_els_names; std::map > neg_moles; std::map els; -double *y, *Ct1, *Ct2, *l_tk_x2, **A, **mixf; +double *Ct2, *l_tk_x2, **A, **LU, **mixf, **mixf_stag; int mixf_comp_size = 0; struct CURRENT_CELLS @@ -35,7 +35,7 @@ struct V_M // For calculating Vinograd and McBain's zero-charge, diffusive tra }; struct CT /* summed parts of V_M and mcd transfer in a timestep for all cells, for free + DL water */ { - LDBLE kgw, dl_s, Dz2c, visc1, visc2, J_ij_sum; + LDBLE kgw, dl_s, Dz2c, Dz2c_stag, visc1, visc2, J_ij_sum; LDBLE A_ij_il, Dz2c_il, mixf_il; int J_ij_count_spec, J_ij_il_count_spec; struct V_M *v_m, *v_m_il; @@ -49,6 +49,7 @@ struct MOLES_ADDED /* total moles added to balance negative conc's */ char *name; LDBLE moles; } *moles_added; + /* ---------------------------------------------------------------------- */ int Phreeqc:: transport(void) @@ -104,9 +105,10 @@ transport(void) sol_D[i].count_exch_spec = 0; sol_D[i].exch_total = 0; sol_D[i].x_max = 0; + sol_D[i].viscos_f = 1.0; + sol_D[i].tk_x = 298.15; sol_D[i].spec = NULL; sol_D[i].spec_size = 0; - sol_D[i].tk_x = 298.15; } //sol_D_dbg = sol_D; @@ -117,9 +119,8 @@ transport(void) { ct[i].kgw = 1.0; ct[i].dl_s = 0.0; - ct[i].Dz2c = 0.0; - ct[i].visc1 = 0.0; - ct[i].visc2 = 0.0; + ct[i].Dz2c = ct[i].Dz2c_stag = 0.0; + ct[i].visc1 = ct[i].visc2 = 0.0; ct[i].J_ij_sum = 0.0; ct[i].A_ij_il = 0.0; ct[i].Dz2c_il = 0.0; @@ -197,6 +198,13 @@ transport(void) use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, k)); if (use.Get_solution_ptr() != NULL) cell_data[k].temp = use.Get_solution_ptr()->Get_tc(); + if (n == 1 && implicit && use.Get_solution_ptr() == NULL) + { + input_error++; + error_string = sformatf( + "Stagnant solution %d not found for implicit diffusion with 1 stagnant layer.\n Please define it, or set -implicit false.", k); + error_msg(error_string, CONTINUE); + } } } @@ -254,14 +262,14 @@ transport(void) { input_error++; error_string = sformatf( - "Interlayer diffusion cannot be calculated implicitly. Please remove -implicit true"); + "Interlayer diffusion cannot be calculated implicitly. Please remove -implicit true."); error_msg(error_string, CONTINUE); } if (implicit && !multi_Dflag) { input_error++; error_string = sformatf( - "Implicit diffusion needs diffusion coefficients for individual species. Please add -multi_d true"); + "Implicit diffusion needs diffusion coefficients for individual species. Please add -multi_d true."); error_msg(error_string, CONTINUE); } if (implicit && !timest) @@ -271,6 +279,12 @@ transport(void) warning_msg(error_string); timest = 1; } + //if (implicit && stag_data->count_stag > 1) + //{ + // error_string = sformatf( + // "Sorry, implicit diffusion can handle only 1 stagnant layer for now. Please remove -implicit true, or set -stagnant 1."); + // //error_msg(error_string, CONTINUE); + //} if (implicit && current_cells == NULL) { current_cells = (struct CURRENT_CELLS *) PHRQ_malloc((size_t) @@ -308,7 +322,7 @@ transport(void) } if (multi_Dflag == TRUE) { - fill_spec(cell_no); + fill_spec(cell_no, 0); } print_punch(i, true); @@ -330,7 +344,7 @@ transport(void) set_and_run_wrapper(k, NOMIX, FALSE, k, 0.0); if (multi_Dflag == TRUE) { - fill_spec(cell_no); + fill_spec(cell_no, i); } print_punch(k, true); saver(); @@ -348,7 +362,7 @@ transport(void) /* multi_D calc's are OK if all cells have the same amount of water */ if (multi_Dflag == TRUE) { - sprintf(token, "multi_D calc's and stagnant: define MIXing factors explicitly, or \n\t adapt the default Dw of -multi_d to match the mobile-immobile exchange factor."); + sprintf(token, "multi_D calc's and stagnant: define MIXing factors explicitly, or \n\t give in -multi_D the Dw used for calculating the mobile-immobile exchange factor."); warning_msg(token); } @@ -541,7 +555,7 @@ transport(void) cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); if (multi_Dflag) - fill_spec(i); + fill_spec(i, i - 1); saver(); } } @@ -553,7 +567,7 @@ transport(void) STOP); } if (implicit) - diffuse_implicit(1, stagkin_time, FALSE); + diffuse_implicit(stagkin_time, stag_data->count_stag); else if (multi_Dflag) multi_D(stagkin_time, 1, FALSE); @@ -586,12 +600,14 @@ transport(void) run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i } else - run_reactions(i, kin_time, DISP, step_fraction); // nsaver = -2 + { + run_reactions(i, kin_time, DISP, step_fraction); // n_saver = -2 + } if (multi_Dflag) - fill_spec(i); + fill_spec(i, 0); /* punch and output file */ - if (ishift == 0 && j == nmix && stag_data->count_stag == 0) + if (ishift == 0 && j == nmix && (stag_data->count_stag == 0 || (implicit && stag_data->count_stag == 1))) print_punch(i, true); if (i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); @@ -740,7 +756,7 @@ transport(void) cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); if (multi_Dflag) - fill_spec(i); + fill_spec(i, i - 1); saver(); } } @@ -762,7 +778,7 @@ transport(void) status(0, token); run_reactions(i, kin_time, NOMIX, step_fraction); if (multi_Dflag == TRUE) - fill_spec(i); + fill_spec(i, i - 1); if (overall_iterations > max_iter) max_iter = overall_iterations; if (nmix == 0 && stag_data->count_stag == 0) @@ -848,7 +864,7 @@ transport(void) cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); if (multi_Dflag) - fill_spec(i); + fill_spec(i, i - 1); saver(); } } @@ -857,9 +873,8 @@ transport(void) if (disp_surf(stagkin_time) == ERROR) error_msg("Error in surface transport, stopping.", STOP); } - if (implicit) - diffuse_implicit(1, stagkin_time, FALSE); + diffuse_implicit(stagkin_time, stag_data->count_stag); else if (multi_Dflag) multi_D(stagkin_time, 1, FALSE); @@ -891,10 +906,12 @@ transport(void) if (i == 0 || i == count_cells + 1) run_reactions(i, kin_time, NOMIX, step_fraction); else + { run_reactions(i, kin_time, DISP, step_fraction); + } if (multi_Dflag == TRUE) - fill_spec(i); - if (j == nmix && stag_data->count_stag == 0) + fill_spec(i, 0); + if (j == nmix && (stag_data->count_stag == 0 || (implicit && stag_data->count_stag == 1))) print_punch(i, true); if (i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); @@ -936,9 +953,11 @@ transport(void) else punch_boolean = FALSE; for (i = 0; i <= count_cells + 1; i++) // allow for stagnant cell mixing with boundary cells + { mix_stag(i, stagkin_time, punch_boolean, step_fraction); + } } - if (j == nmix && stag_data->count_stag > 0) + if (j == nmix && ((stag_data->count_stag > 0/* && !implicit) || (implicit && stag_data->count_stag == 1*/))) { for (n = 1; n <= stag_data->count_stag; n++) { @@ -958,7 +977,7 @@ transport(void) } screen_msg("\n"); - if (implicit) + if (implicit && neg_moles.size()) { std::map > ::iterator it1; std::map ::iterator it2; @@ -1061,23 +1080,36 @@ transport_cleanup(void) } if (implicit) { - y = (LDBLE *)free_check_null(y); - Ct1 = (LDBLE *)free_check_null(Ct1); + int l_stag = (stag_data->count_stag < 2 ? stag_data->count_stag : 0); Ct2 = (LDBLE *)free_check_null(Ct2); l_tk_x2 = (LDBLE *)free_check_null(l_tk_x2); - for (i = 0; i < count_cells + 2; i++) + if (A) { - A[i] = (LDBLE *)free_check_null(A[i]); - mixf[i] = (LDBLE *)free_check_null(mixf[i]); - if (!dV_dcell) + for (i = 0; i < count_cells + 2 + l_stag * count_cells; i++) { - cell_data[i].potV = 0.0; - use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i)); - use.Get_solution_ptr()->Set_potV(0); + A[i] = (LDBLE *)free_check_null(A[i]); + LU[i] = (LDBLE *)free_check_null(LU[i]); + } + } + if (mixf) + { + for (i = 0; i < count_cells + 2; i++) + { + mixf[i] = (LDBLE *)free_check_null(mixf[i]); + if (l_stag) + mixf_stag[i] = (LDBLE *)free_check_null(mixf_stag[i]); + if (!dV_dcell) + { + cell_data[i].potV = 0.0; + use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i)); + use.Get_solution_ptr()->Set_potV(0); + } } } A = (LDBLE **)free_check_null(A); + LU = (LDBLE **)free_check_null(LU); mixf = (LDBLE **)free_check_null(mixf); + mixf_stag = (LDBLE **)free_check_null(mixf_stag); dif_spec_names.clear(); mixf_comp_size = 0; } @@ -1439,7 +1471,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) { if (n == 1) { - if (heat_nmix > 0) + if (heat_nmix > 0 && (!implicit || (implicit && stag_data->count_stag > 1))) { ptr_m = Utilities::Rxn_find(Rxn_solution_map, i); t_imm = @@ -1451,12 +1483,12 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); if (multi_Dflag == TRUE) - fill_spec(cell_no); + fill_spec(cell_no, 0); saver(); cell_no = k; set_and_run_wrapper(k, NOMIX, FALSE, k, 0.0); if (multi_Dflag == TRUE) - fill_spec(cell_no); + fill_spec(cell_no, i); saver(); } /* @@ -1469,18 +1501,17 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) error_msg("Error in surface transport, stopping.", STOP); } - //if (implicit) // work this out appt - // diffuse_implicit(1, stagkin_time, FALSE); - //else if (multi_Dflag) - // multi_D(stagkin_time, 1, FALSE); - if (multi_Dflag == TRUE) - multi_D(1.0, i, TRUE); - set_and_run_wrapper(i, STAG, FALSE, -2, 0.0); - if (multi_Dflag == TRUE) - fill_spec(cell_no); - saver(); // save solution i in -2, original can be used in other stagnant mixes - if (l_punch) - print_punch(i, true); + if (!implicit || (implicit && stag_data->count_stag > 1)) + { + if (multi_Dflag == TRUE) + multi_D(1.0, i, 2); + set_and_run_wrapper(i, STAG, FALSE, -2, 0.0); + if (multi_Dflag == TRUE) + fill_spec(cell_no, 0); + if (l_punch) + print_punch(i, true); + saver(); // save solution i in -2, original can be used in other stagnant mixes + } /* maybe sorb a surface component... */ if (l_punch && change_surf_count) @@ -1501,9 +1532,12 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) } cell_no = k; - run_reactions(k, kin_time, STAG, step_fraction); + if (implicit) + run_reactions(k, kin_time, NOMIX, step_fraction); + else + run_reactions(k, kin_time, STAG, step_fraction); if (multi_Dflag == TRUE) - fill_spec(cell_no); + fill_spec(cell_no, i); saver(); // save solution k in -2 - k, original k can be used in other stagnant mixes /* maybe sorb a surface component... */ @@ -1525,7 +1559,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) done_mixing = true; } - else if (n == 1 && l_punch) + else if (n == 1 && l_punch && !implicit) print_punch(i, false); } @@ -1537,7 +1571,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) if (Utilities::Rxn_find(Rxn_solution_map, k) != 0) { Utilities::Rxn_copy(Rxn_solution_map, -2 - k, k); - if (n == 1) + if (n == 1 && !implicit) Utilities::Rxn_copy(Rxn_solution_map, -2, i); } } @@ -1666,8 +1700,6 @@ init_heat_mix(int l_nmix) { if (implicit) { - // must correct Dw for temp... - fill_spec(i); LDBLE viscos_f; l_heat_nmix = l_nmix; for (i = 1; i <= count_cells + 1; i++) @@ -1706,9 +1738,8 @@ heat_mix(int l_heat_nmix) { for (j = 1; j <= count_cells; j++) temp2[j] = - heat_mix_array[j] * temp1[j - 1] + heat_mix_array[j + 1] * - temp1[j + 1] + (1 - heat_mix_array[j] - - heat_mix_array[j + 1]) * temp1[j]; + heat_mix_array[j] * temp1[j - 1] + heat_mix_array[j + 1] * temp1[j + 1] + + (1 - heat_mix_array[j] - heat_mix_array[j + 1]) * temp1[j]; for (j = 1; j <= count_cells; j++) temp1[j] = temp2[j]; } @@ -1809,7 +1840,7 @@ set_initial_moles(int i) temp_exchange.Set_solution_equilibria(true); temp_exchange.Set_n_solution(i); - cxxExchComp comp; + cxxExchComp comp(this->phrq_io); count_elts = 0; paren_count = 0; strcpy(token, "X"); @@ -1834,7 +1865,7 @@ set_initial_moles(int i) /* ---------------------------------------------------------------------- */ int Phreeqc:: -fill_spec(int l_cell_no) +fill_spec(int l_cell_no, int ref_cell) /* ---------------------------------------------------------------------- */ { /* copy species activities into sol_D.spec... */ @@ -1879,7 +1910,7 @@ fill_spec(int l_cell_no) sol_D[l_cell_no].spec[i].aq_name = NULL; sol_D[l_cell_no].spec[i].type = -1; sol_D[l_cell_no].spec[i].a = 0.0; - sol_D[l_cell_no].spec[i].lm = -3.0; + sol_D[l_cell_no].spec[i].lm = min_dif_LM; sol_D[l_cell_no].spec[i].lg = -0.04; sol_D[l_cell_no].spec[i].c = 0.0; sol_D[l_cell_no].spec[i].z = 0.0; @@ -1889,7 +1920,6 @@ fill_spec(int l_cell_no) sol_D[l_cell_no].count_exch_spec = sol_D[l_cell_no].count_spec = 0; } - sol_D[l_cell_no].tk_x = tk_x; viscos_f = viscos_il_f = 1.0; @@ -1939,19 +1969,21 @@ fill_spec(int l_cell_no) * copy species data */ s_ptr = species_list[i].s; - - if (s_ptr->type == EX && !interlayer_Dflag) - continue; - if (s_ptr->type == SURF) + if (s_ptr->type > HPLUS && + !(s_ptr->type == EX && interlayer_Dflag)) continue; + //if (s_ptr->type == EX && !interlayer_Dflag) + // continue; + //if (s_ptr->type == SURF) + // continue; if (i > 0 && strcmp(s_ptr->name, species_list[i - 1].s->name) == 0) continue; - if (s_ptr == s_h2o) - continue; + //if (s_ptr == s_h2o) + // continue; if (s_ptr->type == EX) { - if (s_ptr->moles > 1e-30) + if (s_ptr->lm > min_dif_LM) { /* find exchanger's name, use only master exchanger 'X' */ if (species_list[i].master_s->secondary != NULL) @@ -2032,9 +2064,9 @@ fill_spec(int l_cell_no) lm = s_ptr->lm; - if (lm > MIN_LM) + if (lm > min_dif_LM) { - if (implicit) + if (implicit/* && l_cell_no < count_cells + 2 */) { name_ret = dif_spec_names.insert(s_ptr->name); loc_spec_names.insert(s_ptr->name); @@ -2066,30 +2098,23 @@ fill_spec(int l_cell_no) i1++; } i3 = i2 - i1; + if (i3 + count_spec + 1 > sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(i3 + count_spec + 1 + size_xt) * sizeof(struct spec)); + if (sol_D[l_cell_no].spec == NULL) + malloc_error(); + sol_D[l_cell_no].spec_size = i3 + count_spec + 1 + size_xt; + } for (i1; i1 < i2; i1++) { - if (i3 > sol_D[l_cell_no].spec_size) - { - sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(i3 + size_xt) * sizeof(struct spec)); - if (sol_D[l_cell_no].spec == NULL) - malloc_error(); - sol_D[l_cell_no].spec_size = i3 + size_xt; - } - memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); + memmove(&sol_D[l_cell_no].spec[i1], &sol_D[ref_cell].spec[i1], sizeof(struct spec)); sol_D[l_cell_no].spec[i1].c = 0.0; sol_D[l_cell_no].spec[i1].a = 0.0; - sol_D[l_cell_no].spec[i1].lm = -3.0; + sol_D[l_cell_no].spec[i1].lm = min_dif_LM; sol_D[l_cell_no].spec[i1].lg = -0.04; loc_spec_names.insert(sol_D[l_cell_no].spec[i1].name); count_spec++; - if (count_spec >= sol_D[l_cell_no].spec_size) - { - sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(count_spec + size_xt) * sizeof(struct spec)); - if (sol_D[l_cell_no].spec == NULL) - malloc_error(); - sol_D[l_cell_no].spec_size = count_spec + size_xt; - } } } } @@ -2130,7 +2155,7 @@ fill_spec(int l_cell_no) diffc_max = sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn); sol_D[l_cell_no].spec[count_spec].erm_ddl = s_ptr->erm_ddl; - if (implicit) + if (implicit/* && l_cell_no < count_cells + 2 */) { if (name_ret.second && l_cell_no) { @@ -2158,7 +2183,7 @@ fill_spec(int l_cell_no) memmove(&sol_D[i1].spec[i2], &sol_D[l_cell_no].spec[i2], sizeof(struct spec)); sol_D[i1].spec[i2].a = 0.0; - sol_D[i1].spec[i2].lm = -3.0; + sol_D[i1].spec[i2].lm = min_dif_LM; sol_D[i1].spec[i2].lg = -0.04; sol_D[i1].spec[i2].c = 0.0; sol_D[i1].count_spec += 1; @@ -2170,7 +2195,7 @@ fill_spec(int l_cell_no) } sol_D[l_cell_no].count_spec = count_spec; sol_D[l_cell_no].count_exch_spec = count_exch_spec; - if (implicit && loc_spec_names.size() < dif_spec_names.size()) + if (implicit/* && l_cell_no < count_cells + 2 */ && loc_spec_names.size() < dif_spec_names.size()) { i3 = (int) dif_spec_names.size(); if (i3 > sol_D[l_cell_no].spec_size) @@ -2182,10 +2207,10 @@ fill_spec(int l_cell_no) } for (i1 = count_spec; i1 < i3; i1++) { - memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); + memmove(&sol_D[l_cell_no].spec[i1], &sol_D[ref_cell].spec[i1], sizeof(struct spec)); sol_D[l_cell_no].spec[i1].c = 0.0; sol_D[l_cell_no].spec[i1].a = 0.0; - sol_D[l_cell_no].spec[i1].lm = -3.0; + sol_D[l_cell_no].spec[i1].lm = min_dif_LM; sol_D[l_cell_no].spec[i1].lg = -0.04; sol_D[l_cell_no].count_spec += 1; @@ -2200,15 +2225,16 @@ fill_spec(int l_cell_no) // error_msg(error_string, CONTINUE); //} } - if (implicit && l_cell_no) + if (implicit && l_cell_no/* && l_cell_no < count_cells + 2 */) { for (i = 0; i < count_spec; i++) { - if (strcmp(sol_D[l_cell_no].spec[i].name, sol_D[l_cell_no - 1].spec[i].name)) + //name = sol_D[l_cell_no].spec[i].name; + if (strcmp(sol_D[l_cell_no].spec[i].name, sol_D[ref_cell].spec[i].name)) { error_string = sformatf( "Implicit diffusion: species %s in cell %d differs from species %s in previous cells.", - sol_D[l_cell_no].spec[i].name, l_cell_no, sol_D[l_cell_no - 1].spec[i].name); + sol_D[l_cell_no].spec[i].name, l_cell_no, sol_D[ref_cell].spec[i].name); error_msg(error_string, CONTINUE); } } @@ -2230,50 +2256,72 @@ sort_species_name(const void *ptr1, const void *ptr2) } /* ---------------------------------------------------------------------- */ void Phreeqc:: -diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) +diffuse_implicit(LDBLE DDt, int stagnant) /* ---------------------------------------------------------------------- */ { // The structure comes from example 11.3 in A&P and Appelo 2017, CCR 101, 102: // Ct2 is implicitly solved for the individual species, corrected for electro-neutrality, this always gives Ct2 > 0. + // Version 3.5.2 includes diffusion in 1 stagnant layer in the implicit calc'n. // Transport of aqueous species is summarized into master species. // With electro-migration, transport of anions and cations is calculated in opposite directions since (sign) J = - z * dV. // Only available moles are transported, thus are > 0, but if concentrations oscillate, // change max_mixf in input file: -implicit true 1 # max_mixf = 1 (default). - int i, ifirst, ilast, cp, comp; - double mfr, mfr1, max_b = 0, b, grad, dVc, j_0e; - cxxSolution *sptr1, *sptr2; + int i, icell, cp, comp; + // ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells); + int ifirst, ilast; + int i_1, i0, i1, i2; + double mfr, mfr1, max_b = 0, b, grad, dVc, j_0e, min_dif_M = pow(10, min_dif_LM); + LDBLE dum1, dum2, dum_stag, min_mol; + + LDBLE dum = 0; + cxxSurfaceCharge * charge_ptr = NULL; + if (stagnant > 1) + stagnant = 0; + cxxSolution *sptr1, *sptr2, *sptr_stag; + std::vector n_solution; + std::vector fraction; + std::map > ::iterator it1; + std::map ::iterator it2; + // c= count_cells, c1=(c+1)= end boundary-cell, c_1=(c-1), c2=(c+2)= first stagnant cell, cc1=(c+c+1)= last stagnant cell + 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; if (heat_nmix) comp += 1; - if (y == NULL) - y = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); - if (y == NULL) malloc_error(); - if (Ct1 == NULL) - Ct1 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); - if (Ct1 == NULL) malloc_error(); if (Ct2 == NULL) - Ct2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + Ct2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2 + stagnant * count_cells) * sizeof(LDBLE)); if (Ct2 == NULL) malloc_error(); if (l_tk_x2 == NULL) - l_tk_x2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + l_tk_x2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2 + stagnant * count_cells) * sizeof(LDBLE)); if (l_tk_x2 == NULL) malloc_error(); - //for (i = 0; i < count_cells + 2; i++) - //{ - // y[i] = 0.0; Ct1[i] = 0.0; Ct2[i] = 0.0; l_tk_x2[i] = 0.0; - //} if (A == NULL) { - A = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2) * sizeof(LDBLE *)); + A = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2 + stagnant * count_cells) * sizeof(LDBLE *)); if (A == NULL) malloc_error(); - for (i = 0; i < count_cells + 2; i++) + for (i = 0; i < count_cells + 2 + stagnant * count_cells; i++) { - A[i] = (LDBLE *)PHRQ_malloc((size_t)3 * sizeof(LDBLE)); + if (stagnant) + A[i] = (LDBLE *)PHRQ_calloc((size_t)(2 * count_cells + 2), sizeof(LDBLE)); + else + A[i] = (LDBLE *)PHRQ_malloc((size_t)3 * sizeof(LDBLE)); if (A[i] == NULL) malloc_error(); } } + if (LU == NULL) + { + LU = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2 + stagnant * count_cells) * sizeof(LDBLE *)); + if (LU == NULL) malloc_error(); + for (i = 0; i < count_cells + 2 + stagnant * count_cells; i++) + { + if (stagnant) + LU[i] = (LDBLE *)PHRQ_calloc((size_t)(2 * count_cells + 2), sizeof(LDBLE)); + else + LU[i] = (LDBLE *)PHRQ_malloc((size_t)3 * sizeof(LDBLE)); + if (LU[i] == NULL) malloc_error(); + } + } if (mixf == NULL) { mixf = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2) * sizeof(LDBLE *)); @@ -2283,6 +2331,18 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) mixf[i] = NULL; } } + if (stagnant) + { + if (mixf_stag == NULL) + { + mixf_stag = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2) * sizeof(LDBLE *)); + if (mixf_stag == NULL) malloc_error(); + for (i = 0; i < count_cells + 2; i++) + { + mixf_stag[i] = NULL; + } + } + } if (comp + 2 > mixf_comp_size) { for (i = 0; i < count_cells + 2; i++) @@ -2292,27 +2352,95 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) else mixf[i] = (LDBLE *)PHRQ_realloc(mixf[i], (size_t)(comp + 2) * sizeof(LDBLE)); if (mixf[i] == NULL) malloc_error(); - mixf_comp_size = comp + 2; + if (stagnant) + { + if (mixf_stag[i] == NULL) + mixf_stag[i] = (LDBLE *)PHRQ_malloc((size_t)(comp + 2) * sizeof(LDBLE)); + else + mixf_stag[i] = (LDBLE *)PHRQ_realloc(mixf_stag[i], (size_t)(comp + 2) * sizeof(LDBLE)); + if (mixf_stag[i] == NULL) malloc_error(); + for (cp = 0; cp < comp; cp++) + mixf_stag[i][cp] = 0; + } } + mixf_comp_size = comp + 2; } - //for (i = 0; i < count_cells + 2; i++) - //{ - // for (int i1 = 0; i1 < 3; i1++) - // A[i][i1] = 0.0; - // for (int i1 = 0; i1 < comp; i1++) - // mixf[i][i1] = 0.0; - //} if (dV_dcell) find_current = 1; // obtain b_ij... - for (i = 0; i <= count_cells; i++) + dummy = default_Dw * pow(multi_Dpor, multi_Dn) * multi_Dpor; + for (i = 0; i <= count_cells + 1; i++) { if (!heat_nmix) - l_tk_x2[i] = (sol_D[i].tk_x + sol_D[i + 1].tk_x) / 2; - b = find_J(i, i + 1, 1, 1, stagnant); - if (b > max_b) - max_b = b; + { + if (i <= count_cells) + l_tk_x2[i] = (sol_D[i].tk_x + sol_D[i + 1].tk_x) / 2; + //if (stagnant && i > 0 && i <= c) + // l_tk_x2[i + c1] = (sol_D[i].tk_x + sol_D[i + c1].tk_x) / 2; + } + if (stagnant) + { + use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i)); + if (use.Get_mix_ptr() == NULL) + { + i1 = 0; + if (i == 0 && (sptr1 = Utilities::Rxn_find(Rxn_solution_map, c2)) != NULL) + { + // the 1st boundary solution diffuses into 1st stagnant cell (=c2) by the ratio immobile_kgw in cell c2 / mobile_kgw in cell 1 + b = find_J(0, 1, 1, 1, 0); + ct[c2].kgw = sptr1->Get_mass_water(); + for (cp = 0; cp < comp; cp++) + { + if (!heat_nmix || cp < comp - 1) + mixf_stag[i][cp] = DDt * ct[i].v_m[cp].b_ij * ct[c2].kgw / ct[1].kgw; + else if (heat_nmix && cp == comp - 1) + mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[i].spec[0].Dwt; + } + i1 = c2; + } + else if (i == c1 && (sptr1 = Utilities::Rxn_find(Rxn_solution_map, cc1)) != NULL) + { + // the end boundary solution diffuses into stagnant cell cc1 + b = find_J(c, c1, 1, 1, 0); + ct[cc1].kgw = sptr1->Get_mass_water(); + for (cp = 0; cp < comp; cp++) + { + if (!heat_nmix || cp < comp - 1) + mixf_stag[i][cp] = DDt * ct[c].v_m[cp].b_ij * ct[cc1].kgw / ct[c].kgw; + else if (heat_nmix && cp == comp - 1) + mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[c].spec[0].Dwt; + } + i1 = cc1; + } + if (i1) + { + b = find_J(i, i1, 1, 1, 1); + } + } + else + { + use.Get_mix_ptr()->Vectorize(n_solution, fraction); + for (i1 = 0; i1 < (int) use.Get_mix_ptr()->Get_mixComps().size(); i1++) + { + if (n_solution[i1] > count_cells && n_solution[i1] <= cc1) + { + mfr = fraction[i1] / nmix; + b = find_J(i, n_solution[i1], mfr, 1, 1); + mfr *= (2.0 / dummy); + for (cp = 0; cp < comp; cp++) + { + if (!heat_nmix || cp < comp - 1) + mixf_stag[i][cp] = ct[i].v_m[cp].b_ij * mfr; + else if (heat_nmix && cp == comp - 1) + mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[i].spec[0].Dwt; + } + } + } + } + } + if (i <= count_cells) + b = find_J(i, i + 1, 1, 1, 0); } if (dV_dcell) find_current = 0; @@ -2322,12 +2450,12 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) for (cp = 0; cp < comp; cp++) { - for (i = 0; i <= count_cells + 1; i++) + for (i = 0; i <= cc1; i++) { if (heat_nmix && cp == comp - 1) - Ct2[i] = Ct1[i] = sol_D[i].tk_x; + Ct2[i] = sol_D[i].tk_x; else - Ct2[i] = Ct1[i] = sol_D[i].spec[cp].c; + Ct2[i] = sol_D[i].spec[cp].c; } // fill coefficient matrix A ... // boundary cells ... @@ -2338,51 +2466,136 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) } else { - mfr = mixf[0][cp] = DDt * ct[0].v_m[cp].b_ij; // maybe add retardation appt + mfr = mixf[0][cp] = DDt * ct[0].v_m[cp].b_ij; mfr1 = mixf[1][cp] = DDt * ct[1].v_m[cp].b_ij; } if (bcon_first == 2) { - A[0][1] = 1; A[0][2] = 0; A[1][0] = 0; A[1][1] = 1 + mfr1; A[1][2] = -mfr1; - } - else - { - if (dV_dcell) + if (stagnant) { - A[0][1] = 1 + mfr; A[0][2] = -mfr; + A[0][0] = A[c2][c2] = 1; + if (mixf_stag[1][cp]) + { + A[1][1] += mixf_stag[1][cp]; + A[1][c2] = A[c2][1] = -mixf_stag[1][cp]; + A[c2][c2] += mixf_stag[1][cp]; + } } else { A[0][1] = 1; A[0][2] = 0; } + } + else + { A[1][0] = -mfr; A[1][1] = 1 + mfr + mfr1; A[1][2] = -mfr1; + if (stagnant) + { + if (dV_dcell) + { + A[0][0] = 1 + mfr; A[0][1] = -mfr; + } + else + A[0][0] = 1; + A[c2][c2] = 1; + if (mixf_stag[0][cp]) + { + if (dV_dcell) + { + A[0][0] += mixf_stag[0][cp]; A[0][c2] = -mixf_stag[0][cp]; + } + A[c2][0] = -mixf_stag[0][cp]; + A[c2][c2] += mixf_stag[0][cp]; + } + if (mixf_stag[1][cp]) + { + A[1][1] += mixf_stag[1][cp]; + A[1][c2] = A[c2][1] = -mixf_stag[1][cp]; + A[c2][c2] += mixf_stag[1][cp]; + } + } + else + { + if (dV_dcell) + { + A[0][1] = 1 + mfr; A[0][2] = -mfr; + } + else + { + A[0][1] = 1; A[0][2] = 0; + } + } } if (heat_nmix && cp == comp - 1) { - mfr = mixf[count_cells - 1][cp] = heat_mix_array[count_cells - 1]; + mfr = mixf[c_1][cp] = heat_mix_array[c_1]; mfr1 = mixf[count_cells][cp] = heat_mix_array[count_cells]; } else { - mfr = mixf[count_cells - 1][cp] = DDt * ct[count_cells - 1].v_m[cp].b_ij; + mfr = mixf[c_1][cp] = DDt * ct[c_1].v_m[cp].b_ij; mfr1 = mixf[count_cells][cp] = DDt * ct[count_cells].v_m[cp].b_ij; } if (bcon_last == 2) { - A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0; - A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; - } - else - { - A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1; - if (dV_dcell) + if (stagnant) { - A[count_cells + 1][0] = -mfr1; A[count_cells + 1][1] = 1 + mfr1; + A[count_cells][c_1] = -mfr; A[count_cells][count_cells] = 1 + mfr; + A[c1][c1] = A[cc1][cc1] = 1; + if (mixf_stag[count_cells][cp]) + { + A[count_cells][count_cells] += mixf_stag[count_cells][cp]; + A[count_cells][cc1] = A[cc1][count_cells] = -mixf_stag[count_cells][cp]; + A[cc1][cc1] += mixf_stag[count_cells][cp]; + } } else { - A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0; + A[c1][0] = 0; A[c1][1] = 1; + } + } + else + { + if (stagnant) + { + A[count_cells][c_1] = -mfr; A[count_cells][count_cells] = 1 + mfr + mfr1; + A[count_cells][c1] = -mfr1; + if (dV_dcell) + { + A[c1][count_cells] = -mfr1; A[c1][c1] = 1 + mfr1; + } + else + A[c1][c1] = 1; + A[cc1][cc1] = 1; + if (mixf_stag[count_cells][cp]) + { + A[count_cells][count_cells] += mixf_stag[count_cells][cp]; + A[count_cells][cc1] = A[cc1][count_cells] = -mixf_stag[count_cells][cp]; + A[cc1][cc1] += mixf_stag[count_cells][cp]; + } + if (mixf_stag[c1][cp]) + { + if (dV_dcell) + { + A[c1][c1] += mixf_stag[c1][cp]; A[c1][cc1] = -mixf_stag[c1][cp]; + } + A[cc1][c1] = -mixf_stag[c1][cp]; + A[cc1][cc1] += mixf_stag[c1][cp]; + } + } + else + { + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1; + if (dV_dcell) + { + A[c1][0] = -mfr1; A[c1][1] = 1 + mfr1; + } + else + { + A[c1][0] = 0; A[c1][1] = 1; + } } } // inner cells ... @@ -2398,25 +2611,139 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) mfr = mixf[i - 1][cp] = DDt * ct[i - 1].v_m[cp].b_ij; mfr1 = mixf[i][cp] = DDt * ct[i].v_m[cp].b_ij; } - A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1; + if (stagnant && mixf_stag[i][cp]) + { + A[i][i - 1] = -mfr; A[i][i] = 1 + mfr + mfr1 + mixf_stag[i][cp]; A[i][i + 1] = -mfr1; + A[i + c1][i + c1] = 1 + mixf_stag[i][cp]; + A[i][i + c1] = A[i + c1][i] = -mixf_stag[i][cp]; + } + else + { + A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1; + } } - // decompose A in LU : store L in A[..][0..1] and U in A[..][2] ... - for (i = 1; i <= count_cells + 1; i++) + if (stagnant) + // decompose full A in LU, bypass A[][] = 0... { - A[i - 1][2] = A[i - 1][2] / A[i - 1][1]; - A[i][1] = A[i][1] - A[i][0] * A[i - 1][2]; + LDBLE s; + LU[0][0] = A[0][0]; LU[1][0] = A[1][0]; LU[c2][0] = A[c2][0]; + LU[0][1] = A[0][1] / LU[0][0]; LU[0][c2] = A[0][c2] / LU[0][0]; + for (i = 1; i <= cc1; i++) + { + i_1 = i - 1; + if (i < c2) + { + LU[i][i] = A[i][i] - LU[i][i_1] * LU[i_1][i]; + LU[i + 1][i] = A[i + 1][i]; + LU[i][i + 1] = A[i][i + 1] / LU[i][i]; + if (i == 1) + { + LU[c2][1] = A[c2][1] - LU[c2][0] * LU[0][1]; + LU[1][c2] = (A[1][c2] - LU[1][0] * LU[0][c2]) / LU[1][1]; + } + else + { + for (i0 = 0; i0 <= i_1; i0++) + { + if (i0 == i_1 && i0 < c) + { + LU[c2 + i0][i] = A[c2 + i0][i]; + LU[i][c2 + i0] = A[i][c2 + i0] / LU[i][i]; + } + else if (i0 < c) + { + LU[c2 + i0][i] = -LU[c2 + i0][i_1] * LU[i_1][i]; + LU[i][c2 + i0] = -LU[i][i_1] * LU[i_1][c2 + i0] / LU[i][i]; + } + if (i > c && i0 == c_1) + { + LU[c2 + i0][i] += A[c2 + i0][i]; + LU[i][c2 + i0] += A[i][c2 + i0] / LU[i][i]; + } + } + } + } + else + { + // the i = count_cells + 2 and higher rows are filled from i0 onwards, need to subtract the L*U terms for the ith column of L. + for (i0 = i; i0 <= cc1; i0++) + { + s = 0; + if (i0 == i) + { + i2 = (i == c2 ? 0 : i - c1); + for (i1 = i2; i1 <= i_1; i1++) + s += LU[i0][i1] * LU[i1][i]; + LU[i][i] = A[i][i] - s; + i2 += (i2 == 0 ? 2 : 1); + } + else + { + for (i1 = i2; i1 <= i_1; i1++) + s += LU[i0][i1] * LU[i1][i]; + LU[i0][i] = -s; + } + } + // and the ith row of U... + if (i == cc1) + break; + for (i0 = i + 1; i0 <= cc1; i0++) + { + s = 0; + if (i0 == i + 1) + i2 = i - c; + for (i1 = i2; i1 <= i_1; i1++) + s += LU[i][i1] * LU[i1][i0]; + LU[i][i0] = -s / LU[i][i]; + i2++; + } + } + } + Ct2[0] /= LU[0][0]; + for (i = 1; i <= cc1; i++) + { + i_1 = i - 1; + if (i < c2) + s = LU[i][i_1] * Ct2[i_1]; + else + { + i2 = 0; + s = 0; + for (i1 = i2; i1 <= i_1; i1++) + { + s += LU[i][i1] * Ct2[i1]; + i2 += (i2 == 0 ? 2 : 1); + } + } + Ct2[i] = (Ct2[i] - s) / LU[i][i]; + } + for (i = cc; i >= 0; i--) + { + s = 0; + for (i1 = i + 1; i1 <= cc1; i1++) + s += LU[i][i1] * Ct2[i1]; + Ct2[i] -= s; + } } - // solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1 ... - // First, find y in L.y = Ct1 - y[0] = Ct2[0]; - for (i = 1; i <= count_cells + 1; i++) - y[i] = (Ct2[i] - A[i][0] * y[i - 1]) / A[i][1]; - // Now obtain Ct2 in U.Ct2 = y ... - Ct2[count_cells + 1] = y[count_cells + 1]; - for (i = count_cells; i >= 0; i--) - Ct2[i] = y[i] - A[i][2] * Ct2[i + 1]; - - // Moles transported by concentration gradient from cell [i] to [i + 1] go in tot1. + else + { + // decompose A in LU : store L in A[..][0..1] and U in A[..][2] ... + for (i = 1; i <= count_cells + 1; i++) + { + A[i - 1][2] = A[i - 1][2] / A[i - 1][1]; + A[i][1] -= A[i][0] * A[i - 1][2]; + } + // solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1, Ct1 was put in Ct2 ... + // First, find y in L.y = Ct1, put y in Ct2 + Ct2[0] /= A[0][1]; + for (i = 1; i <= count_cells + 1; i++) + Ct2[i] = (Ct2[i] - A[i][0] * Ct2[i - 1]) / A[i][1]; + // Now obtain Ct2 in U.Ct2 = 'y' ... + for (i = count_cells; i >= 0; i--) + Ct2[i] -= A[i][2] * Ct2[i + 1]; + } + // Moles transported by concentration gradient from cell [i] to [i + 1] go in tot1, + // moles by stagnant exchange from cell [i] to [i + c1] go in tot_stag // Correct for electro-neutrality... for (i = ifirst; i <= ilast + 1; i++) { @@ -2426,28 +2753,63 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) cell_data[i].temp = Ct2[i] - 273.15; sptr1 = Utilities::Rxn_find(Rxn_solution_map, i); sptr1->Set_tc(Ct2[i] - 273.15); + if (stagnant && i > 0 && i < c1 && mixf_stag[i][cp]) + { + i1 = i + c1; + cell_data[i1].temp = Ct2[i1] - 273.15; + sptr2 = Utilities::Rxn_find(Rxn_solution_map, i1); + sptr2->Set_tc(Ct2[i1] - 273.15); + } continue; } - else if (i == ilast + 1) + else if (i == ilast + 1 && (!stagnant || !mixf_stag[i][cp])) continue; - if (cp == 0) - current_cells[i].ele = current_cells[i].dif = 0; - grad = Ct2[i + 1] - Ct2[i]; - ct[i].J_ij[cp].name = sol_D[i].spec[cp].name; - ct[i].J_ij[cp].tot1 = -mixf[i][cp] * grad; - ct[i].J_ij[cp].charge = ct[i].v_m[cp].z; - if (ct[i].v_m[cp].z) + if (!cp) { - if (i == 0) - ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i + 1]; - else if (i == ilast) - ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i]; + ct[i].J_ij_sum = 0; + ct[i].Dz2c_stag = 0; + if (i < ilast + 1) + current_cells[i].ele = current_cells[i].dif = 0; + } + for (i0 = stagnant; i0 >= 0; i0--) + { + i1 = (i0 == 0 ? i + 1 : i == 0 ? c2 : i == c1 ? cc1 : i + c1); + if (Utilities::Rxn_find(Rxn_solution_map, i1) == NULL) + continue; + ct[i].J_ij[cp].name = sol_D[i].spec[cp].name; + 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; else - ct[i].v_m[cp].zc = ct[i].v_m[cp].z * (Ct2[i] + Ct2[i + 1]) / 2; - - current_cells[i].dif += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot1; - current_cells[i].ele -= ct[i].v_m[cp].z * mixf[i][cp] * ct[i].v_m[cp].zc; + ct[i].J_ij[cp].tot_stag = -mixf_stag[i][cp] * grad; + if (ct[i].v_m[cp].z) + { + if (i == 0) + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i1]; + else if (i == ilast) + { + if (!i0) + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i]; + else + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i1]; + } + else + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * (Ct2[i] + Ct2[i1]) / 2; + if (i0) + { + mixf_stag[i][cp] *= ct[i].v_m[cp].zc; + ct[i].J_ij_sum += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot_stag; + ct[i].Dz2c_stag -= ct[i].v_m[cp].z * mixf_stag[i][cp]; + } + else if (i < ilast + 1) + { + mixf[i][cp] *= ct[i].v_m[cp].zc; + current_cells[i].dif += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot1; + current_cells[i].ele -= ct[i].v_m[cp].z * mixf[i][cp]; + } + } } } // define element list @@ -2499,7 +2861,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) else { current_x = 1e-10 / F_C_MOL; - cell_data[ifirst].potV = 1e-8; + cell_data[ifirst].potV = 0e-8; j_0e = current_x - current_cells[ifirst].dif; } dVc = j_0e * current_cells[ifirst].R; @@ -2521,28 +2883,31 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) for (cp = 0; cp < comp; cp++) { - for (i = ifirst; i <= ilast; i++) + if (!ct[ifirst].v_m[cp].z) + continue; + for (i = ifirst; i <= ilast + stagnant; i++) { - ct[i].J_ij_sum = 0; - if (!ct[i].v_m[cp].z) - continue; - 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] * ct[i].v_m[cp].zc * dVc; - // With dVc known, charge transport = 0 if dV_dcell = 0... - if (!dV_dcell) - ct[i].J_ij_sum += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot1; + if (i < ilast + 1) + { + 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; + } + if (stagnant && ct[i].Dz2c_stag) + { + ct[i].J_ij[cp].tot_stag += (mixf_stag[i][cp] * ct[i].J_ij_sum / ct[i].Dz2c_stag); + } } } current_A = current_x / DDt * F_C_MOL; - for (i = ifirst; i <= ilast + 1; i++) + for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++) { - // preserve the potentials... - sptr1 = Utilities::Rxn_find(Rxn_solution_map, i); - sptr1->Set_potV(cell_data[i].potV); - if (i == ilast + 1) - continue; - + if (i <= ilast + 1) + { + // preserve the potentials... + sptr1 = Utilities::Rxn_find(Rxn_solution_map, i); + sptr1->Set_potV(cell_data[i].potV); + } // Translate transport of the solute species into master species... ct[i].count_m_s = count_m_s; if (ct[i].m_s_size == 0 && ct[i].m_s != NULL) @@ -2560,23 +2925,20 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) if (ct[i].m_s == NULL) malloc_error(); std::set ::iterator it = dif_els_names.begin(); - for (int i1 = 0; i1 < count_m_s; i1++) + for (i1 = 0; i1 < count_m_s; i1++) { - ct[i].m_s[i1].tot1 = 0.0; + ct[i].m_s[i1].tot1 = ct[i].m_s[i1].tot2 = ct[i].m_s[i1].tot_stag = 0.0; ct[i].m_s[i1].charge = 0.0; ct[i].m_s[i1].name = (*it).c_str(); it++; } - fill_m_s(ct[i].J_ij, ct[i].J_ij_count_spec, i, 0); + fill_m_s(ct[i].J_ij, ct[i].J_ij_count_spec, i, stagnant); } /* - * 3. find the solutions, add or subtract the moles... + * 3. find the solutions in the column, add or subtract the moles... */ - std::map > ::iterator it1; - std::map ::iterator it2; cxxNameDouble::iterator kit; - int icell, if1, il1, incr, length, length2; - LDBLE dum1, dum2, min_mol; + int if1, il1, incr; for (cp = 0; cp < count_m_s; cp++) { if (!dV_dcell || dV_dcell * ct[1].m_s[cp].charge <= 0) @@ -2589,35 +2951,76 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) } for (icell = if1; icell != il1; icell += incr) { - min_mol = 1e-12 * ct[icell].kgw; - if (fabs(ct[icell].m_s[cp].tot1) < min_mol) - continue; + min_mol = min_dif_M * ct[icell].kgw; dum1 = dum2 = 0; sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell); sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1); + if (stagnant) + { + i1 = (icell == 0 ? c1 + 1 : icell == c1 ? cc1 : icell + c1); + sptr_stag = Utilities::Rxn_find(Rxn_solution_map, i1); + } + else + sptr_stag = NULL; + if (!strcmp(ct[icell].m_s[cp].name, "H")) { + dummy = ct[icell].m_s[cp].tot1; if (dV_dcell || (icell > 0 && icell <= ilast)) - sptr1->Set_total_h(sptr1->Get_total_h() - ct[icell].m_s[cp].tot1); + sptr1->Set_total_h(sptr1->Get_total_h() - dummy); if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) - sptr2->Set_total_h(sptr2->Get_total_h() + ct[icell].m_s[cp].tot1); + sptr2->Set_total_h(sptr2->Get_total_h() + dummy); + if (sptr_stag) + { + dummy = ct[icell].m_s[cp].tot_stag; + if (dV_dcell || (icell > 0 && icell <= ilast)) + sptr1->Set_total_h(sptr1->Get_total_h() - dummy); + if (icell == c) + { + // mix the boundary solution with the stagnant column end-cell + dummy += ct[icell + 1].m_s[cp].tot_stag; + if (dV_dcell) + sptr2->Set_total_h(sptr2->Get_total_h() - ct[icell + 1].m_s[cp].tot_stag); + } + sptr_stag->Set_total_h(sptr_stag->Get_total_h() + dummy); + } continue; } if (!strcmp(ct[icell].m_s[cp].name, "O")) { + dummy = ct[icell].m_s[cp].tot1; if (dV_dcell || (icell > 0 && icell <= ilast)) - sptr1->Set_total_o(sptr1->Get_total_o() - ct[icell].m_s[cp].tot1); + sptr1->Set_total_o(sptr1->Get_total_o() - dummy); if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) - sptr2->Set_total_o(sptr2->Get_total_o() + ct[icell].m_s[cp].tot1); + sptr2->Set_total_o(sptr2->Get_total_o() + dummy); + if (sptr_stag) + { + dummy = ct[icell].m_s[cp].tot_stag; + if (dV_dcell || (icell > 0 && icell <= ilast)) + sptr1->Set_total_o(sptr1->Get_total_o() - dummy); + if (icell == c) + { + dummy += ct[icell + 1].m_s[cp].tot_stag; + if (dV_dcell) + sptr2->Set_total_o(sptr2->Get_total_o() - ct[icell + 1].m_s[cp].tot_stag); + } + 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). + //{ + // sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum); + // sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum); + // if (sptr_stag) + // sptr_stag->Set_cb(sptr_stag->Get_cb() + ct[i1].J_ij_sum); + //} continue; } - // see if icell - incr has negative moles, subtract it, so that it is not transported... + // 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; - // the negative moles are canceled... neg_moles.erase(it1); els.erase(it2); neg_moles.insert(std::make_pair(icell - incr, els)); @@ -2625,34 +3028,27 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name]; if (!dum1) { - /* Get moles from redox states... */ - length = (int)strlen(ct[icell].m_s[cp].name); - for (kit = sptr1->Get_totals().begin(); kit != sptr1->Get_totals().end(); kit++) - { - length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); - if (length == length2 && !strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) - { - dum1 += kit->second; kit->second = 0; - } - } + 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) { - length = (int)strlen(ct[icell].m_s[cp].name); - for (kit = sptr2->Get_totals().begin(); kit != sptr2->Get_totals().end(); kit++) - { - length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); - if (length == length2 && !strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) - { - dum2 += kit->second; kit->second = 0; - } - } + 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 && @@ -2665,31 +3061,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell); if (s_ptr) { - cxxSurfaceCharge * charge_ptr = NULL; - for (size_t j = 0; j < s_ptr->Get_surface_charges().size(); j++) - { - if (s_ptr->Get_dl_type() == cxxSurface::DONNAN_DL) - { - charge_ptr = &(s_ptr->Get_surface_charges()[j]); - for (kit = charge_ptr->Get_diffuse_layer_totals().begin(); kit != charge_ptr->Get_diffuse_layer_totals().end(); kit++) - { - //if (strcmp(kit->first.c_str(), "H") == 0 || strcmp(kit->first.c_str(), "O") == 0) - // continue; - if (strcmp(kit->first.c_str(), ct[icell].m_s[cp].name) == 0) - { - dummy = ct[icell].m_s[cp].tot1 - dum1 + min_mol; - if (kit->second > dummy) - { - dum1 += dummy; kit->second -= dummy; - } - else - { - dum1 += kit->second; kit->second = 0; - } - } - } - } - } + dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 - dum1 + min_mol); } } sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1; @@ -2706,31 +3078,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1); if (s_ptr) { - cxxSurfaceCharge * charge_ptr = NULL; - for (size_t j = 0; j < s_ptr->Get_surface_charges().size(); j++) - { - if (s_ptr->Get_dl_type() == cxxSurface::DONNAN_DL) - { - charge_ptr = &(s_ptr->Get_surface_charges()[j]); - for (kit = charge_ptr->Get_diffuse_layer_totals().begin(); kit != charge_ptr->Get_diffuse_layer_totals().end(); kit++) - { - //if (strcmp(kit->first.c_str(), "H") == 0 || strcmp(kit->first.c_str(), "O") == 0) - // continue; - if (strcmp(kit->first.c_str(), ct[icell].m_s[cp].name) == 0) - { - dummy = ct[icell].m_s[cp].tot2 + dum2 - min_mol; - if (kit->second > -dummy) - { - dum2 -= dummy; kit->second += dummy; - } - else - { - dum2 += kit->second; kit->second = 0; - } - } - } - } - } + 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) @@ -2740,32 +3088,214 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) 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) - { - els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1)); - neg_moles.insert(std::make_pair(icell, els)); + if (dum1 < -min_mol && incr > 0) + { + 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; } 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) + { + 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; } - if (dum2 < -min_mol && incr < 0) + + if (sptr_stag) { - els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2)); - neg_moles.insert(std::make_pair(icell + 1, els)); + 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; + if (icell == c) + { + 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) + { + 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)); + } + 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) + { + 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; } + //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); + // if (sptr_stag) + // sptr_stag->Set_cb(sptr_stag->Get_cb() + ct[i1].J_ij_sum); + //} } } return; } +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +moles_from_redox_states(cxxSolution *sptr, const char *name) +/* ---------------------------------------------------------------------- */ +{ + int length, length2; + cxxNameDouble::iterator kit; + LDBLE dum = 0; + + length = (int)strlen(name); + for (kit = sptr->Get_totals().begin(); kit != sptr->Get_totals().end(); kit++) + { + length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); + if (length == length2 && !strncmp(name, kit->first.c_str(), length2)) + { + dum += kit->second; kit->second = 0; + } + } + return(dum); +} + +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +add_MCD_moles(LDBLE moles, LDBLE min_mol, int i, cxxSolution *sptr, const char *name) +/* ---------------------------------------------------------------------- */ +{ + cxxNameDouble::iterator kit; + std::map > ::iterator it1; + std::map ::iterator it2; + LDBLE dum = sptr->Get_totals()[name]; + + if (!dum) + dum = moles_from_redox_states(sptr, name); + if ((it1 = neg_moles.find(i)) != neg_moles.end() + && (it2 = (els = it1->second).find(name)) != els.end()) + { + dum += it2->second; + neg_moles.erase(it1); + els.erase(it2); + neg_moles.insert(std::make_pair(i, els)); + } + dum += moles; + if (dum < -min_mol) + { + if (ct[i].dl_s) + { + cxxSurface *su_ptr = Utilities::Rxn_find(Rxn_surface_map, i); + if (su_ptr != NULL) + dum += moles_from_donnan_layer(su_ptr, name, -dum + min_mol); + } + } + sptr->Get_totals()[name] = (dum > 0 ? dum : 0e-16); + if (dum < -min_mol) + { + els.insert(std::make_pair(name, dum)); + if ((it1 = neg_moles.find(i)) != neg_moles.end()) + neg_moles.erase(it1); + neg_moles.insert(std::make_pair(i, els)); + } + return(dum); +} + +/* ---------------------------------------------------------------------- */ +LDBLE Phreeqc:: +moles_from_donnan_layer(cxxSurface *sptr, const char *name, LDBLE moles_needed) +/* ---------------------------------------------------------------------- */ +{ + cxxNameDouble::iterator kit; + LDBLE dum = 0; + cxxSurfaceCharge * charge_ptr = NULL; + + for (size_t j = 0; j < sptr->Get_surface_charges().size(); j++) + { + if (sptr->Get_dl_type() == cxxSurface::DONNAN_DL) + { + charge_ptr = &(sptr->Get_surface_charges()[j]); + for (kit = charge_ptr->Get_diffuse_layer_totals().begin(); kit != charge_ptr->Get_diffuse_layer_totals().end(); kit++) + { + if (strcmp(kit->first.c_str(), "H") == 0 || strcmp(kit->first.c_str(), "O") == 0) + continue; + if (strcmp(kit->first.c_str(), name) == 0) + { + if (kit->second > moles_needed) + { + dum += moles_needed; kit->second -= moles_needed; + } + else + { + dum += kit->second; kit->second = 0; + } + } + } + } + } + return(dum); +} /* ---------------------------------------------------------------------- */ int Phreeqc:: @@ -3014,8 +3544,8 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) // first_c = mobile_cell; // allow for stagnant cell mixing with boundary cell 0 for (i = first_c; i <= last_c2; i++) { - if (stagnant && i > first_c && i <= count_cells + first_c) - continue; + //if (stagnant && i > first_c && i <= count_cells + first_c) + // continue; use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i)); if (!use.Get_solution_ptr()) @@ -3137,8 +3667,6 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant) for (j = 0; j < l_J_ij_count_spec; j++) { - if (fabs(l_J_ij[j].tot1) < 1e-15) - continue; { char * temp_name = string_duplicate(l_J_ij[j].name); ptr = temp_name; @@ -3146,7 +3674,7 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant) get_elts_in_species(&ptr, 1); free_check_null(temp_name); } - if (implicit && !stagnant) + if (implicit && stagnant < 2) { for (k = 0; k < count_elts; k++) { @@ -3160,10 +3688,10 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant) else fraction = 1; ct[icell].m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; - //if (!strcmp(elt_list[k].elt->name, "H") || !strcmp(elt_list[k].elt->name, "O")) - // break; ct[icell].m_s[l].charge *= (1 - fraction); ct[icell].m_s[l].charge += fraction * l_J_ij[j].charge; + if (stagnant) + ct[icell].m_s[l].tot_stag += elt_list[k].coef * l_J_ij[j].tot_stag; break; } } @@ -3272,13 +3800,13 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) std::vector::iterator it_sc; std::vector s_com_p; - ct[icell].dl_s = dl_aq1 = dl_aq2 = 0.0; + ct[icell].dl_s = ct[jcell].dl_s = dl_aq1 = dl_aq2 = 0.0; if (dV_dcell && !find_current) goto dV_dcell2; /* check for immediate return and interlayer diffusion calcs... */ - ct[icell].J_ij_sum = 0.0; + //ct[icell].J_ij_sum = 0.0; ct[icell].J_ij_count_spec = 0; if (!il_calcs) { @@ -3393,7 +3921,7 @@ 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 = dl_aq2 / t_aq2; + ct[icell].dl_s = ct[jcell].dl_s = dl_aq2 / t_aq2; if (il_calcs) { @@ -3566,9 +4094,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) for (i = 0; i < k; i++) { - ct[icell].J_ij[i].tot1 = 0.0; + ct[icell].J_ij[i].tot1 = ct[icell].J_ij[i].tot_stag = 0.0; ct[icell].v_m[i].grad = 0.0; -// ct[icell].v_m[i].D = 0.0; + ct[icell].v_m[i].D = 1.0; // used for gamma correction ct[icell].v_m[i].z = 0.0; ct[icell].v_m[i].c = 0.0; ct[icell].v_m[i].zc = 0.0; @@ -3922,7 +4450,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) 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) + 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) @@ -3932,9 +4460,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) } else { - if (icell == 0) + 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) + else if (icell == count_cells && jcell == count_cells + 1) ct[icell].v_m[k].b_ij = b_i; } @@ -3945,18 +4473,18 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) //if (fabs(ddlm) > 1e-10) // ct[icell].v_m[k].grad *= (1 + (sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg) / ddlm); ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; - dum = sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg; - if (fabs(dum * sol_D[jcell].spec[j].lg) > 5e-3 && fabs(ddlm) > 0.025) + dum1 = sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg; + dum = 1; + if (ddlm) + dum += dum1 / ddlm; + if (dum > 0.2 && dum < 2.25) { - dum = (1 + dum / ddlm); ct[icell].v_m[k].grad *= dum; if (implicit) - ct[icell].v_m[k].b_ij *= dum; + ct[icell].v_m[k].D = dum; } - - if (ct[icell].v_m[k].b_ij > max_b) - max_b = ct[icell].v_m[k].b_ij; - + //if (ct[icell].v_m[k].b_ij > max_b) + // max_b = ct[icell].v_m[k].b_ij; k++; } if (i < i_max) @@ -3965,20 +4493,20 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) j++; } } - if (implicit && !stagnant) - return(max_b); - /* - * fill in J_ij... - */ if (!dV_dcell && !ct[icell].Dz2c) k = 0; ct[icell].J_ij_count_spec = i_max = k; ct[icell].J_ij_il_count_spec = k_il; + if (icell == count_cells) + ct[jcell].J_ij_count_spec = ct[icell].J_ij_count_spec; + if (implicit && stagnant < 2) + return(max_b); + /* + * fill in J_ij... + */ if (dV_dcell) { - //if (transport_step >= 100) // debug... - // icell = icell; current_cells[icell].ele = current_cells[icell].dif = 0; for (i = 0; i < ct[icell].J_ij_count_spec; i++) { @@ -3999,7 +4527,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) // voltage was adapted to give equal current in the cells. dV_dcell2: - ct[icell].J_ij_sum = 0; + //ct[icell].J_ij_sum = 0; Sum_zM = 0.0; for (i = 0; i < ct[icell].J_ij_count_spec; i++) @@ -4012,21 +4540,23 @@ dV_dcell2: ct[icell].J_ij[i].tot1 = -ct[icell].v_m[i].grad; 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; + } if (stagnant) ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * 2 * mixf; else ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * DDt; 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; + //ct[icell].J_ij_sum += ct[icell].v_m[i].z * ct[icell].J_ij[i].tot1; } - // assure that icell has dl water when checking negative conc's in MCD + // 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; if (dV_dcell) { - ct[icell].J_ij_sum = 0; + //ct[icell].J_ij_sum = 0; dV = cell_data[jcell].potV - cell_data[icell].potV; dum = dV * F_Re3 / tk_x2; // perhaps adapt dV for getting equal current... @@ -4051,15 +4581,10 @@ dV_dcell2: ct[icell].J_ij[i].tot1 -= ct[icell].v_m[i].b_ij * ct[icell].v_m[i].zc * dum * DDt; 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; + //ct[icell].J_ij_sum += ct[icell].v_m[i].z * ct[icell].J_ij[i].tot1; } current_A = current_x * F_C_MOL; } - //if (fabs(ct[icell].J_ij_sum / DDt - current_x) > 1e-13 && (dV_dcell || (icell > 0 || jcell >= count_cells))) // appt - //{ - // sprintf(token, "Moles charge added in cell %d is: %.3e.\n", icell, ct[icell].J_ij_sum / DDt - current_x); - // screen_msg(token); - //} /* * calculate interlayer mass transfer... */ @@ -4079,7 +4604,7 @@ dV_dcell2: ct[icell].J_ij_il[i].tot1 *= ct[icell].mixf_il; else ct[icell].J_ij_il[i].tot1 *= ct[icell].A_ij_il * DDt; - ct[icell].J_ij_sum += ct[icell].v_m_il[i].z * ct[icell].J_ij_il[i].tot1; + //ct[icell].J_ij_sum += ct[icell].v_m_il[i].z * ct[icell].J_ij_il[i].tot1; ct[icell].J_ij_il[i].tot2 = ct[icell].J_ij_il[i].tot1; // ct[icell].J_ij_il[i].charge = ct[icell].v_m_il[i].z; // appt not used now }