Merge commit '63108ac76d6fa7c378c21690b6342333c4cb446c'

src/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 <<

src/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);
 	}

src/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;

src/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;

src/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);
+	void diffuse_implicit(LDBLE DDt, int stagnant);
-	int fill_spec(int cell_no);
+	int fill_spec(int cell_no, int ref_cell);
-	void define_ct_structures(void);
+	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 */

src/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);
 		}

src/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;

src/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;

src/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();

src/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);
 		}

src/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;

src/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.2");
-	len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1");
 #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

src/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
 				}

src/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

src/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);
 	}

src/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

src/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",
+				"\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());
+				(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
 			error_msg(error_string, STOP);
 		}
 	}

src/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.",
+					"Bad RK steps > %d in cell %d. Please decrease (time)step or increase -bad_step_max.",
-					kinetics_ptr->Get_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();

src/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);

src/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 *) malloc(sizeof(PHRQMemHeader) + size * num);
-	//p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
+	p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
 
 	if (p == NULL)
 		return NULL;

src/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
  */

src/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);
+		cxxPPassemblage * pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, use.Get_n_pp_assemblage_user());
-		cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr;
+		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
  */

src/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
 				*/

src/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)
 	{