debugging and optimizing

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

Exchange.cxx

@@ -37,6 +37,7 @@ cxxNumKeyword(io)
 {
 	this->n_user = this->n_user_end = l_n_user;
 	this->pitzer_exchange_gammas = true;
+	this->new_def = false;
 //
 //   Mix exchangers
 //

NameDouble.cxx

@@ -297,6 +297,7 @@ cxxNameDouble::Simplify_redox(void) const
 	new_totals.type = cxxNameDouble::ND_ELT_MOLES;
 	{
 		std::string current_ename;
+		std::string const *ename_ptr;
 		cxxNameDouble::const_iterator it;
 
 		// make list of elements in new_totals
@@ -305,23 +306,24 @@ cxxNameDouble::Simplify_redox(void) const
 			current_ename = it->first;
 			if (it->first.size() < 4)
 			{
-				current_ename = it->first;
+				ename_ptr = &(it->first);
 			}
 			else
 			{
-				indexCh = current_ename.find("(");
+				indexCh = it->first.find("(");
 				if (indexCh != std::string::npos)
 				{
-					current_ename = current_ename.substr(0, indexCh);
+					current_ename = it->first.substr(0, indexCh);
+					ename_ptr = &(current_ename);
 				}
 				else
 				{
-					current_ename = it->first;
+					ename_ptr = &(it->first);
 				}
 			}
 			if (current_ename == "H" || current_ename == "O" || current_ename == "Charge")
 				continue;
-			new_totals[current_ename] = 0;
+			new_totals[*ename_ptr] = 0;
 		}
 	}
 
@@ -330,25 +332,27 @@ cxxNameDouble::Simplify_redox(void) const
 		cxxNameDouble::const_iterator old_it = nd.begin();
 		cxxNameDouble::iterator new_it = new_totals.begin();
 		std::string old_ename;
+		std::string const *old_ename_ptr;
 		while (old_it != nd.end() && new_it != new_totals.end())
 		{
 			if (old_it->first.size() < 4)
 			{
-				old_ename = old_it->first;
+				old_ename_ptr = &old_it->first;
 			}
 			else
 			{
 				indexCh = old_it->first.find("(");
 				if (indexCh != std::string::npos)
 				{
-					old_ename = old_ename.substr(0, indexCh);
+					old_ename = old_it->first.substr(0, indexCh);
+					old_ename_ptr = &old_ename;
 				}
 				else
 				{
-					old_ename = old_it->first;
+					old_ename_ptr = &old_it->first;
 				}
 			}
-			int j = strcmp(new_it->first.c_str(), old_ename.c_str());
+			int j = strcmp(new_it->first.c_str(), old_ename_ptr->c_str());
 			if (j < 0)
 			{
 				new_it++;
@@ -407,6 +411,35 @@ cxxNameDouble::Multiply_activities_redox(std::string str, LDBLE f)
 	std::string redox_name = str;
 	redox_name.append("(");
 
+	for (it = this->begin(); it != this->end(); it++)
+	{
+		if (str[0] > it->first[0]) continue;
+		if (it->first == str)
+		{
+			// Found exact match
+			it->second += lg_f;
+		}
+		else 
+		{
+			// no exact match, current is element name, need to find all valences
+			if (strstr(it->first.c_str(), redox_name.c_str()) == it->first.c_str())
+			{
+				it->second += lg_f;
+			}
+		}
+		if (str[0] < it->first[0]) break;
+	}
+}
+#ifdef SKIP
+void 
+cxxNameDouble::Multiply_activities_redox(std::string str, LDBLE f)
+{
+	// update original master_activities using just computed factors
+	cxxNameDouble::iterator it;
+	LDBLE lg_f = log10(f);
+	std::string redox_name = str;
+	redox_name.append("(");
+
 	for (it = this->begin(); it != this->end(); it++)
 	{
 		if (it->first == str)
@@ -424,6 +457,7 @@ cxxNameDouble::Multiply_activities_redox(std::string str, LDBLE f)
 		}
 	}
 }
+#endif
 LDBLE
 cxxNameDouble::Get_total_element(const char *string) const
 {
@@ -539,7 +573,6 @@ cxxNameDouble::merge_redox(const cxxNameDouble & source)
 			// Put in elt name
 			(*this)[elt_name] = sit->second;
 		}
-
 	}
 }
 struct DblCmp {     

Solution.cxx

@@ -288,6 +288,350 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
 
 	return;
 }
+
+#ifdef USE_REVISED_READ_RAW
+void
+cxxSolution::read_raw(CParser & parser, bool check)
+{
+
+	// Used if it is modify
+	cxxNameDouble original_totals = this->totals;
+	cxxNameDouble original_activities(this->master_activity);
+
+	this->master_activity.clear();
+
+	std::istream::pos_type ptr;
+	std::istream::pos_type next_char;
+	std::string token;
+	int opt_save;
+
+	// Read solution number and description
+	this->read_number_description(parser.line());
+
+	opt_save = CParser::OPT_ERROR;
+	bool tc_defined(false);
+	bool ph_defined(false);
+	bool pe_defined(false);
+	bool mu_defined(false);
+	bool ah2o_defined(false);
+	bool total_h_defined(false);
+	bool total_o_defined(false);
+	bool cb_defined(false);
+	bool mass_water_defined(false);
+	bool total_alkalinity_defined(false);
+
+	for (;;)
+	{
+		int opt = parser.get_option(vopts, next_char);
+		if (opt == CParser::OPT_DEFAULT)
+		{
+			opt = opt_save;
+		}
+
+		switch (opt)
+		{
+		case CParser::OPT_EOF:
+			break;
+		case CParser::OPT_KEYWORD:
+			break;
+		case CParser::OPT_DEFAULT:
+		case CParser::OPT_ERROR:
+			opt = CParser::OPT_EOF;
+			parser.error_msg("Unknown input in SOLUTION_RAW keyword.",
+							 PHRQ_io::OT_CONTINUE);
+			parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE);
+			continue;
+
+		case 0:				// totals
+			{
+				cxxNameDouble temp_totals;
+				if (temp_totals.read_raw(parser, next_char) !=	CParser::PARSER_OK)
+				{
+					parser.incr_input_error();
+					parser.
+						error_msg("Expected element name and moles for totals.",
+						PHRQ_io::OT_CONTINUE);
+				}
+				else
+				{
+					this->totals.merge_redox(temp_totals);
+				}
+			}
+			opt_save = 0;
+			break;
+
+		case 1:				// activities
+			if (this->master_activity.read_raw(parser, next_char) !=
+				CParser::PARSER_OK)
+			{
+				parser.incr_input_error();
+				parser.
+					error_msg
+					("Expected species name and log activity for activities.",
+					 PHRQ_io::OT_CONTINUE);
+			}
+			opt_save = 1;
+			break;
+
+		case 2:				// gammas
+			if (this->species_gamma.read_raw(parser, next_char) !=
+				CParser::PARSER_OK)
+			{
+				parser.incr_input_error();
+				parser.
+					error_msg
+					("Expected species name and activity coefficient for gammas.",
+					 PHRQ_io::OT_CONTINUE);
+			}
+			opt_save = 2;
+			break;
+
+		case 3:				// isotope
+			{
+				std::string name;
+				if (!(parser.get_iss() >> name))
+				{
+					parser.incr_input_error();
+					parser.error_msg("Expected character value for isotope name.",
+									 PHRQ_io::OT_CONTINUE);
+				}
+				else
+				{
+					cxxSolutionIsotope iso(this->Get_io());
+					iso.Set_isotope_name(name.c_str());
+					iso.read_raw(parser, check);
+					this->isotopes[name] = iso;
+				}
+			}
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 4:				// temp
+		case 5:				// tc_avoid_conflict_with_technetium
+		case 6:				// temperature                  
+			if (!(parser.get_iss() >> this->tc))
+			{
+				this->tc = 25.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for temperature.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			tc_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 7:				// ph
+			if (!(parser.get_iss() >> this->ph))
+			{
+				this->ph = 7.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for pH.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			ph_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 8:				// pe
+			if (!(parser.get_iss() >> this->pe))
+			{
+				this->pe = 4.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for pe.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			pe_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 9:				// mu
+		case 10:				// ionic_strength
+			if (!(parser.get_iss() >> this->mu))
+			{
+				this->mu = 1e-7;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for ionic strength.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			mu_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 11:				// ah2o
+		case 12:				// activity_water
+			if (!(parser.get_iss() >> this->ah2o))
+			{
+				this->ah2o = 1.0;
+				parser.incr_input_error();
+				parser.
+					error_msg("Expected numeric value for activity of water.",
+							  PHRQ_io::OT_CONTINUE);
+			}
+			ah2o_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 13:				// total_h
+			if (!(parser.get_iss() >> this->total_h))
+			{
+				this->total_h = 111.1;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for total hydrogen.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			total_h_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 14:				// total_o
+			if (!(parser.get_iss() >> this->total_o))
+			{
+				this->total_o = 55.55;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for total oxygen.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			total_o_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 15:				// mass_water
+		case 16:				// mass_h2o
+			if (!(parser.get_iss() >> this->mass_water))
+			{
+				this->mass_water = 1.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for mass of water.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			mass_water_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 17:				// total_alkalinity
+		case 18:				// total_alk
+			if (!(parser.get_iss() >> this->total_alkalinity))
+			{
+				this->total_alkalinity = 0;
+				parser.incr_input_error();
+				parser.
+					error_msg("Expected numeric value for total_alkalinity.",
+							  PHRQ_io::OT_CONTINUE);
+			}
+			total_alkalinity_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+
+		case 19:				// cb
+		case 20:				// charge_balance
+			if (!(parser.get_iss() >> this->cb))
+			{
+				this->cb = 0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for charge balance.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			cb_defined = true;
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+		case 21:				// density
+			if (!(parser.get_iss() >> this->density))
+			{
+				this->density = 1.0;
+				parser.incr_input_error();
+				parser.error_msg("Expected numeric value for density.",
+								 PHRQ_io::OT_CONTINUE);
+			}
+			opt_save = CParser::OPT_DEFAULT;
+			break;
+		}
+
+		if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
+			break;
+	}
+	if (check)
+	{
+		// all members must be defined
+		if (tc_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Temp not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (ph_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("pH not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (pe_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("pe not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (mu_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Ionic strength not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (ah2o_defined == false)
+		{
+			parser.incr_input_error();
+			parser.
+				error_msg("Activity of water not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (total_h_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Total hydrogen not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (total_o_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Total oxygen not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (cb_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Charge balance not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (mass_water_defined == false)
+		{
+			parser.incr_input_error();
+			parser.error_msg("Temp not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+		if (total_alkalinity_defined == false)
+		{
+			parser.incr_input_error();
+			parser.
+				error_msg("Total alkalinity not defined for SOLUTION_RAW input.",
+				PHRQ_io::OT_CONTINUE);
+		}
+	}
+
+	// Update activities
+	
+	if (original_activities.size() > 0)
+	{
+		cxxNameDouble new_activities = this->master_activity;
+		this->master_activity = original_activities;
+		this->Update_activities(original_totals);
+		cxxNameDouble::iterator it = new_activities.begin();
+		for ( ; it != new_activities.end(); it++)
+		{
+			this->master_activity[it->first] = it->second;
+		}
+	}
+	return;
+}
+#endif
 void
 cxxSolution::read_raw(CParser & parser, bool check)
 {
@@ -641,6 +985,7 @@ cxxSolution::read_raw(CParser & parser, bool check)
 
 	return;
 }
+
 void
 cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble &const_nd)
 {
@@ -652,12 +997,13 @@ cxxSolution::Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble
 	this->Update(const_nd);
 }
 void
-cxxSolution::Update(const cxxNameDouble &const_nd)
+cxxSolution::Update_activities(const cxxNameDouble &original_tot)
 {
 	// Totals and activities of solution are updated
 	// nd does not have H, O, charge
-	cxxNameDouble simple_original = this->totals.Simplify_redox();
-	cxxNameDouble simple_new = const_nd.Simplify_redox();
+	cxxNameDouble simple_original = original_tot.Simplify_redox();
+	// totals after read
+	cxxNameDouble simple_new = this->totals.Simplify_redox();
 
 	// make factors 
 	cxxNameDouble factors;
@@ -725,15 +1071,84 @@ cxxSolution::Update(const cxxNameDouble &const_nd)
 			}
 		}
 	}
+}
+void
+cxxSolution::Update(const cxxNameDouble &const_nd)
+{
+	// const_nd is a list of new totals, assumed to be inclusive of all elements
+	// Totals and activities of solution are updated
+	// nd does not have H, O, charge
+	cxxNameDouble simple_original = this->totals.Simplify_redox();
+	cxxNameDouble simple_new = const_nd.Simplify_redox();
 
-	// update other totals
+	cxxNameDouble factors;
 	{
-		cxxNameDouble::iterator it;
-		for (it = simple_new.begin(); it  != simple_new.end(); it++)
+		// make factors 
+		cxxNameDouble::iterator it = simple_new.begin();
+		cxxNameDouble::iterator jit = simple_original.begin();
+		while (it != simple_new.end() && jit != simple_original.end())
 		{
-			simple_original[it->first] = it->second;
+			int j = strcmp(it->first.c_str(), jit->first.c_str());
+			if (j < 0)
+			{
+				it++;
+			}
+			else if (j == 0)
+			{
+				if (jit->second != it->second)
+				{
+					if (it->second > 0 && jit->second > 0)
+					{
+						factors[it->first] = log10(jit->second / it->second);
+					}
+				}
+				it++;
+				jit++;
+			}
+			else
+			{
+				jit++;
+			}
+		}
+
+		// simple_new now has factors for master activities
+		// Now add log factors to log activities
+		{
+			cxxNameDouble::iterator activity_it = this->master_activity.begin();
+			cxxNameDouble::iterator factors_it = factors.begin();
+			std::string activity_ename;
+			std::basic_string < char >::size_type indexCh;
+			while (activity_it != master_activity.end() && factors_it != factors.end())
+			{
+				activity_ename = activity_it->first;
+				if (activity_ename.size() > 3)
+				{
+					indexCh = activity_ename.find("(");
+					if (indexCh != std::string::npos)
+					{
+						activity_ename = activity_ename.substr(0, indexCh);
+					}
+				}
+				int j = strcmp(factors_it->first.c_str(), activity_ename.c_str());
+				if (j < 0)
+				{
+					factors_it++;
+				}
+				else if (j == 0)
+				{
+					activity_it->second += factors_it->second;
+					activity_it++;
+				}
+				else 
+				{
+					activity_it++;
+				}
+			}
 		}
 	}
+
+	// update totals
+	this->totals = simple_new;
 }
 #ifdef SKIP
 void

Solution.h

@@ -104,7 +104,7 @@ class cxxSolution:public cxxNumKeyword
 	//void Simplify_totals();
 	void Update(const cxxNameDouble &nd);
 	void Update(LDBLE h_tot, LDBLE o_tot, LDBLE charge, const cxxNameDouble &nd);
-	void Update_activities(const cxxNameDouble &old_tot, const cxxNameDouble &new_tot);
+	void Update_activities(const cxxNameDouble &original_tot);
 
 #ifdef USE_MPI
 	void mpi_pack(std::vector < int >&ints, std::vector < LDBLE >&doubles);

phreeqc/step.cpp

@@ -308,7 +308,7 @@ xsolution_zero(void)
 	ph_x = 0.0;
 	solution_pe_x = 0.0;
 	mu_x = 0.0;
-	ah2o_x = 1.0;
+	ah2o_x = 0.0;
 	density_x = 0.0;
 	total_h_x = 0.0;
 	total_o_x = 0.0;