Bugs with kinetics in PHAST.

The formula was not correct when converted to a phreeqcpp class, only the last instance of an element was kept and the rest were lost. Modified NameDouble to accumulate the sum of coefficients for each element when converting from a PHREEQC structure.

The order of the kinetics components was not maintained when converted to a phreeqcpp class. Components were in alphabetical order. This posed a problem with get and put statements because the gets could end up before the puts. Changed from map of components to list of components, which should maintain the order correctly.

Problems with modify solution and redox elements. Laurin pointed out that adding N may leave all the N(x) in place, which increases the total N. Revised read_solution to merge the valence states. If N defined, then all N(x) are removed. If N() defined, then N is removed. Still need to test.


git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@3927 1feff8c3-07ed-0310-ac33-dd36852eb9cd

KineticsComp.cxx

@@ -59,7 +59,7 @@ cxxKineticsComp::~cxxKineticsComp()
 }
 
 struct kinetics_comp *
-	cxxKineticsComp::cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::map < std::string, cxxKineticsComp >
+	cxxKineticsComp::cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::list < cxxKineticsComp >
 											   &el)
 		//
 		// Builds kinetics_comp structure from of cxxKineticsComp 
@@ -72,20 +72,20 @@ struct kinetics_comp *
 		P_INSTANCE_POINTER malloc_error();
 
 	int i = 0;
-	for (std::map < std::string, cxxKineticsComp >::iterator it = el.begin();
+	for (std::list < cxxKineticsComp >::iterator it = el.begin();
 		 it != el.end(); ++it)
 	{
-		if ((*it).second.rate_name.size() == 0)
+		if ((*it).rate_name.size() == 0)
 			kinetics_comp_ptr[i].rate_name = NULL;
 		else
-			kinetics_comp_ptr[i].rate_name = P_INSTANCE_POINTER string_hsave((*it).second.rate_name.c_str());
-		kinetics_comp_ptr[i].list = (*it).second.namecoef.name_coef(P_INSTANCE);
-		kinetics_comp_ptr[i].count_list = (int) (*it).second.namecoef.size();
-		kinetics_comp_ptr[i].tol = (*it).second.tol;
-		kinetics_comp_ptr[i].m = (*it).second.m;
+			kinetics_comp_ptr[i].rate_name = P_INSTANCE_POINTER string_hsave((*it).rate_name.c_str());
+		kinetics_comp_ptr[i].list = (*it).namecoef.name_coef(P_INSTANCE);
+		kinetics_comp_ptr[i].count_list = (int) (*it).namecoef.size();
+		kinetics_comp_ptr[i].tol = (*it).tol;
+		kinetics_comp_ptr[i].m = (*it).m;
 		kinetics_comp_ptr[i].initial_moles = 0.;
-		kinetics_comp_ptr[i].m0 = (*it).second.m0;
-		kinetics_comp_ptr[i].moles = (*it).second.moles;
+		kinetics_comp_ptr[i].m0 = (*it).m0;
+		kinetics_comp_ptr[i].moles = (*it).moles;
 		kinetics_comp_ptr[i].count_c_params = 0;
 		kinetics_comp_ptr[i].c_params = NULL;
 /*
@@ -93,15 +93,15 @@ struct kinetics_comp *
                 kinetics_comp_ptr[i].d_params               =  NULL;
 */
 
-		kinetics_comp_ptr[i].count_d_params = (int) (*it).second.d_params.size();
+		kinetics_comp_ptr[i].count_d_params = (int) (*it).d_params.size();
 		kinetics_comp_ptr[i].d_params = NULL;
-		if ((*it).second.d_params.size() > 0)
+		if ((*it).d_params.size() > 0)
 		{
 			kinetics_comp_ptr[i].d_params =	(double *)
-				P_INSTANCE_POINTER PHRQ_malloc((size_t) ((*it).second.d_params.size() * sizeof(double)));
+				P_INSTANCE_POINTER PHRQ_malloc((size_t) ((*it).d_params.size() * sizeof(double)));
 			if (kinetics_comp_ptr[i].d_params == NULL)
 				P_INSTANCE_POINTER malloc_error();
-			std::copy((*it).second.d_params.begin(), (*it).second.d_params.end(),
+			std::copy((*it).d_params.begin(), (*it).d_params.end(),
 					  kinetics_comp_ptr[i].d_params);
 		}
 		i++;

KineticsComp.h

@@ -17,7 +17,7 @@ public:
 	cxxKineticsComp(struct kinetics_comp *);
 	~cxxKineticsComp();
 
-	static struct kinetics_comp *cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::map < std::string, cxxKineticsComp > &el);
+	static struct kinetics_comp *cxxKineticsComp2kinetics_comp(PHREEQC_PTR_ARG_COMMA std::list < cxxKineticsComp > &el);
 
 	void dump_xml(std::ostream & os, unsigned int indent = 0) const;
 

NameDouble.cxx

@@ -114,6 +114,7 @@ cxxNameDouble::cxxNameDouble(struct master_activity *ma, int count,
 	}
 	this->type = type;
 }
+/*
 cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count)
 		//
 		// constructor for cxxNameDouble from list of elt_list
@@ -128,6 +129,30 @@ cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count)
 	}
 	this->type = ND_NAME_COEF;
 }
+*/
+cxxNameDouble::cxxNameDouble(struct name_coef *nc, int count)
+		//
+		// constructor for cxxNameDouble from list of elt_list
+		//
+{
+	int i;
+	for (i = 0; i < count; i++)
+	{
+		if (nc[i].name == NULL)
+			continue;
+
+		if ((*this).find(nc[i].name) == (*this).end())
+		{
+			(*this)[nc[i].name] = nc[i].coef;
+		}
+		else
+		{
+			(*this)[nc[i].name] = (*this).find(nc[i].name)->second + nc[i].coef;
+		}
+
+	}
+	this->type = ND_NAME_COEF;
+}
 
 cxxNameDouble::~cxxNameDouble()
 {
@@ -462,7 +487,70 @@ cxxNameDouble::multiply(double extensive)
 		it->second *= extensive;
 	}
 }
+void
+cxxNameDouble::merge_redox(const cxxNameDouble & source)
+//
+// Merges source into this
+// Accounts for possible conflicts between redox state and
+// totals
+//
+{
+	for (cxxNameDouble::const_iterator sit = source.begin(); sit != source.end(); sit++)
+	{
+		
+		std::string redox_name = sit->first;
+		std::string elt_name;
+		size_t pos = redox_name.find("(");
 
+		bool redox;
+		if (pos != std::string::npos) 
+		{
+			redox = true;
+			elt_name = redox_name.substr(0, pos - 1);
+		}
+		else
+		{
+			redox = false;
+			elt_name = redox_name;
+		}
+		if (redox)
+		{
+			// Remove elt_name, if present
+			if ((*this).find(elt_name) != (*this).end())
+			{
+				(*this).erase((*this).find(elt_name));
+			}
+			// Put in redox name
+			(*this)[redox_name] = sit->second;
+		}
+		else
+		{
+			std::string substring;
+			substring.append(elt_name);
+			substring.append("(");
+
+			// Remove all redox
+			bool deleted = true;
+			while (deleted)
+			{
+				deleted = false;
+				cxxNameDouble::iterator current = (*this).begin();
+				for ( ; current != (*this).end(); current++)
+				{
+					if (current->first.find(substring) == 0)
+					{
+						(*this).erase(current);
+						deleted = true;
+						break;
+					}
+				}
+			}
+			// Put in elt name
+			(*this)[elt_name] = sit->second;
+		}
+
+	}
+}
 #ifdef USE_MPI
 void
 cxxNameDouble::mpi_pack(std::vector < int >&ints,

NameDouble.h

@@ -66,6 +66,9 @@ class cxxNameDouble:public
 	void
 	multiply(double factor);
 
+	void
+	merge_redox(const cxxNameDouble & source);
+
 	void
 	insert(char *str, double d)
 	{

Solution.cxx

@@ -630,13 +630,27 @@ cxxSolution::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check)
 			continue;
 
 		case 0:				// totals
-			if (this->totals.read_raw(P_INSTANCE_COMMA parser, next_char) !=
-				CParser::PARSER_OK)
+			//if (this->totals.read_raw(P_INSTANCE_COMMA parser, next_char) !=
+			//	CParser::PARSER_OK)
+			//{
+			//	parser.incr_input_error();
+			//	parser.
+			//		error_msg("Expected element name and moles for totals.",
+			//				  CParser::OT_CONTINUE);
+			//}
 			{
-				parser.incr_input_error();
-				parser.
-					error_msg("Expected element name and moles for totals.",
-							  CParser::OT_CONTINUE);
+				cxxNameDouble temp_totals;
+				if (temp_totals.read_raw(P_INSTANCE_COMMA parser, next_char) !=	CParser::PARSER_OK)
+				{
+					parser.incr_input_error();
+					parser.
+						error_msg("Expected element name and moles for totals.",
+						CParser::OT_CONTINUE);
+				}
+				else
+				{
+					this->totals.merge_redox(temp_totals);
+				}
 			}
 			opt_save = 0;
 			break;

cxxKinetics.cxx

@@ -66,7 +66,8 @@ totals(kinetics_ptr->totals)
 	{
 		cxxKineticsComp ec(&(kinetics_ptr->comps[i]));
 		std::string str(ec.get_rate_name());
-		this->kineticsComps[str] = ec;
+		//this->kineticsComps[str] = ec;
+		this->kineticsComps.push_back(ec);
 	}
 
 	// steps
@@ -231,12 +232,12 @@ cxxKinetics::dump_raw(std::ostream & s_oss, unsigned int indent) const
 	s_oss << "-cvode_order       " << this->cvode_order << std::endl;
 
 	// kineticsComps structures
-	for (std::map < std::string, cxxKineticsComp >::const_iterator it =
+	for (std::list < cxxKineticsComp >::const_iterator it =
 		 kineticsComps.begin(); it != kineticsComps.end(); ++it)
 	{
 		s_oss << indent1;
 		s_oss << "-component" << std::endl;
-		(*it).second.dump_raw(s_oss, indent + 2);
+		(*it).dump_raw(s_oss, indent + 2);
 	}
 
 	// totals
@@ -421,6 +422,28 @@ cxxKinetics::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check)
 				reread.set_echo_file(CParser::EO_NONE);
 				reread.set_echo_stream(CParser::EO_NONE);
 
+				std::list < cxxKineticsComp >::iterator kit;
+				bool found = false;
+				for (kit = this->kineticsComps.begin(); kit != this->kineticsComps.end(); kit++)
+				{
+					if (kit->get_rate_name() == ec.get_rate_name())
+					{
+						found = true;
+						break;
+					}
+				}
+				if (found)
+				{
+					kit->read_raw(P_INSTANCE_COMMA reread, false);
+				}
+				else
+				{
+					cxxKineticsComp ec1;
+					ec1.read_raw(P_INSTANCE_COMMA reread, false);
+					std::string str(ec1.get_rate_name());
+					this->kineticsComps.push_back(ec1);
+				}
+				/*
 				if (this->kineticsComps.find(ec.get_rate_name()) != this->kineticsComps.end())
 				{
 					cxxKineticsComp & comp = this->kineticsComps.find(ec.get_rate_name())->second;
@@ -433,6 +456,7 @@ cxxKinetics::read_raw(PHREEQC_PTR_ARG_COMMA CParser & parser, bool check)
 					std::string str(ec1.get_rate_name());
 					this->kineticsComps[str] = ec1;
 				}
+				*/
 			}
 			useLastLine = true;
 			break;
@@ -554,11 +578,18 @@ cxxKinetics::mpi_pack(std::vector < int >&ints,
 {
 	ints.push_back(this->n_user);
 	ints.push_back((int) this->kineticsComps.size());
+	std::list < cxxKineticsComp >::iterator it;
+	for (it = this->kineticsComps.begin(); it != this->kineticsComps.end(); it++)
+	{
+		(*it).mpi_pack(ints, doubles);
+	}
+	/*
 	for (std::map < std::string, cxxKineticsComp >::iterator it =
 		 this->kineticsComps.begin(); it != this->kineticsComps.end(); it++)
 	{
 		(*it).second.mpi_pack(ints, doubles);
 	}
+	*/
 	ints.push_back((int) this->steps.size());
 	for (std::vector < double >::iterator it = this->steps.begin();
 		 it != this->steps.end(); it++)
@@ -588,8 +619,9 @@ cxxKinetics::mpi_unpack(int *ints, int *ii, double *doubles, int *dd)
 	{
 		cxxKineticsComp kc;
 		kc.mpi_unpack(ints, &i, doubles, &d);
-		std::string str(kc.get_rate_name());
-		this->kineticsComps[str] = kc;
+		//std::string str(kc.get_rate_name());
+		//this->kineticsComps[str] = kc;
+		this->kineticsComps.push_back(kc);
 	}
 	n = ints[i++];
 	this->steps.clear();
@@ -664,21 +696,31 @@ cxxKinetics::add(const cxxKinetics & addee, double extensive)
 	if (extensive == 0.0)
 		return;
 	//std::map < std::string, cxxKineticsComp> kineticsComps;
-	for (std::map < std::string, cxxKineticsComp >::const_iterator itadd =
+	for (std::list < cxxKineticsComp >::const_iterator itadd =
 		 addee.kineticsComps.begin(); itadd != addee.kineticsComps.end();
 		 ++itadd)
 	{
-		std::map < std::string, cxxKineticsComp >::iterator it = this->kineticsComps.find((*itadd).first);
-		if (it != this->kineticsComps.end())
+		bool found(false);
+		std::list <  cxxKineticsComp >::iterator it;
+		for (it = this->kineticsComps.begin(); it != this->kineticsComps.end(); it++)
 		{
-			(*it).second.add((*itadd).second, extensive);
+			if ((*it).get_rate_name() == (*itadd).get_rate_name())
+			{
+				found = true;
+				break;
+			}
+		}
+
+		//std::map < std::string, cxxKineticsComp >::iterator it = this->kineticsComps.find((*itadd).first);
+		if (found)
+		{
+			(*it).add((*itadd), extensive);
 		}
 		else
 		{
-			cxxKineticsComp entity = (*itadd).second;
+			cxxKineticsComp entity = (*itadd);
 			entity.multiply(extensive);
-			std::string str(entity.get_rate_name());
-			this->kineticsComps[str] = entity;
+			this->kineticsComps.push_back(entity);
 		}
 	}
 	//std::vector<double> steps;

cxxKinetics.h

@@ -43,7 +43,8 @@ class cxxKinetics:public cxxNumKeyword
 	void add(const cxxKinetics & addee, double extensive);
 
   protected:
-	std::map < std::string, cxxKineticsComp > kineticsComps;
+	//std::map < std::string, cxxKineticsComp > kineticsComps;
+	std::list < cxxKineticsComp > kineticsComps;
 	std::vector < double >steps;
 	cxxNameDouble totals;
 	double step_divide;