Updated with phreeqc files.

Main difference is the logic for rebuilding a model. Should
run faster now, but has slight differences in example files.

Imported changes from ph2orch.

Plan to make phreeqcpp a subdirectory of ph2orch/src.




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

ISolution.cxx

@@ -6,13 +6,16 @@
 #endif
 
 #include "ISolution.h"
+#include "ISolutionComp.h"
 #define EXTERNAL extern
 #include "global.h"
 #include "phqalloc.h"
 #include "phrqproto.h"
+#include "output.h"
 
 #include <cassert>     // assert
 #include <algorithm>   // std::sort 
+#include <sstream>
 
 //////////////////////////////////////////////////////////////////////
 // Construction/Destruction
@@ -38,7 +41,8 @@ cxxISolution::cxxISolution(struct solution *solution_ptr)
         // totals
         for (int i = 0; solution_ptr->totals[i].description != NULL; i++) {
                 cxxISolutionComp c(&(solution_ptr->totals[i]));
-                comps.insert(c);
+                //comps.insert(solution_ptr->totals[i].description, c);
+		comps[solution_ptr->totals[i].description] = c;
         }
         default_pe  = solution_ptr->default_pe;
         // pe_data
@@ -68,6 +72,88 @@ struct solution *cxxISolution::cxxISolution2solution()
         soln_ptr->totals = cxxISolutionComp::cxxISolutionComp2conc(this->comps);
         return(soln_ptr);
 }
+void cxxISolution::ConvertUnits()
+  //
+  // Converts from input units to moles per kilogram water
+  //
+{
+  double sum_solutes = 0;
+  // foreach conc
+  std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::iterator iter = this->comps.begin();
+  for(; iter != this->comps.end(); ++iter)
+  {
+    struct master *master_ptr = master_bsearch (iter->first);
+    if (master_ptr != NULL && (master_ptr->minor_isotope == TRUE)) continue;
+    if (iter->second.get_description() == "H(1)" || iter->second.get_description() == "E") continue;
+    if (iter->second.get_input_conc() <= 0.0) continue;
+/*
+*   Convert liters to kg solution
+*/
+    double moles = iter->second.get_input_conc();
+    if (this->units.find("/l") != std::string::npos )
+    {
+	moles /= this->density;
+    }
+/*
+* Convert to moles
+*/
+    //set gfw for element
+    iter->second.set_gfw();
+    // convert to moles
+    if (iter->second.get_units().find("g/") != std::string::npos)
+    {
+      if (iter->second.get_gfw() != 0)
+      {
+	moles /= iter->second.get_gfw();
+      }
+      else 
+      {
+	std::ostringstream oss;
+	oss << "Could not find gfw, " <<  iter->second.get_description();
+	error_msg(oss.str().c_str(), CONTINUE);
+	input_error++;
+      }
+    }
+/*
+*   Convert milli or micro
+*/
+    char c = iter->second.get_units().c_str()[0];
+    if (c == 'm')
+    {
+	moles *= 1e-3;
+    }
+    else if (c == 'u')
+    {
+	moles *= 1e-6;
+    }
+    iter->second.set_moles(moles);
+/*
+*   Sum grams of solute, convert from moles necessary
+*/
+    sum_solutes += moles * (iter->second.get_gfw());
+  }
+/*
+ *   Convert /kgs to /kgw
+ */
+  double mass_water;
+  if ((this->units.find("kgs") != std::string::npos) ||
+    (this->units.find("/l") != std::string::npos))
+  {
+    mass_water = 1.0 - 1e-3 * sum_solutes;
+    for(; iter != this->comps.end(); ++iter)
+    {
+      iter->second.set_moles(iter->second.get_moles() / mass_water);
+    }
+  }
+/*
+ *   Scale by mass of water in solution
+ */
+  mass_water = this->mass_water;
+  for(; iter != this->comps.end(); ++iter)
+  {
+    iter->second.set_moles(iter->second.get_moles() * mass_water);
+  }
+}
 #ifdef SKIP
 cxxISolution& cxxISolution::read(CParser& parser)
 {
@@ -124,8 +210,8 @@ cxxISolution& cxxISolution::read(CParser& parser)
                         break;
                 case CParser::OPTION_ERROR:
                         opt = CParser::OPTION_EOF;
-                        parser.error_msg("Unknown input in SOLUTION keyword.", CParser::OT_CONTINUE);
+                        CParser::error_msg("Unknown input in SOLUTION keyword.", CParser::OT_CONTINUE);
-                        parser.error_msg(parser.line().c_str(), CParser::OT_CONTINUE);
+                        CParser::error_msg(parser.line().c_str(), CParser::OT_CONTINUE);
                         break;
 
                 case 0: // temp

ISolution.h

@@ -39,11 +39,13 @@ public:
 
         //void dump_xml(std::ostream& os, unsigned int indent = 0)const;
 
+	void ConvertUnits();
+
 protected:
         friend class cxxISolutionComp; // for this->pe access
         double density;
         std::string units;
-        std::set<cxxISolutionComp> comps; 
+        std::map<char *, cxxISolutionComp, CHARSTAR_LESS> comps; 
         struct pe_data *pes;
         int default_pe;
 

ISolutionComp.cxx

@@ -10,6 +10,7 @@
 #include "global.h"
 #include "phrqproto.h"
 #include "phqalloc.h"
+#include "output.h"
 
 cxxISolutionComp::cxxISolutionComp(void)
 : description(NULL)
@@ -71,7 +72,7 @@ struct conc *cxxISolutionComp::concarray(std::map <char *, double, CHARSTAR_LESS
         return(c);
 }
 */
-struct conc *cxxISolutionComp::cxxISolutionComp2conc(const std::set <cxxISolutionComp> &totals)
+struct conc *cxxISolutionComp::cxxISolutionComp2conc(const std::map <char *, cxxISolutionComp, CHARSTAR_LESS> &totals)
         // for ISolutions
         // takes a std::vector cxxISolutionComp structures
         // returns list of conc structures
@@ -80,16 +81,16 @@ struct conc *cxxISolutionComp::cxxISolutionComp2conc(const std::set <cxxISolutio
         c = (struct conc *) PHRQ_malloc((size_t) ((totals.size() + 1) * sizeof(struct conc)));
         if (c == NULL) malloc_error();
         int i = 0;
-        for (std::set<cxxISolutionComp>::const_iterator it = totals.begin(); it != totals.end(); ++it) {
+        for (std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::const_iterator it = totals.begin(); it != totals.end(); ++it) {
-                c[i].description         = it->description;
+                c[i].description         = it->second.description;
-                c[i].moles               = it->moles;
+                c[i].moles               = it->second.moles;
-                c[i].input_conc          = it->input_conc;
+                c[i].input_conc          = it->second.input_conc;
-                c[i].units               = it->units;
+                c[i].units               = it->second.units;
-                c[i].equation_name       = it->equation_name;
+                c[i].equation_name       = it->second.equation_name;
-                c[i].phase_si            = it->phase_si;
+                c[i].phase_si            = it->second.phase_si;
-                c[i].n_pe                = it->n_pe;
+                c[i].n_pe                = it->second.n_pe;
-                c[i].as                  = it->as;
+                c[i].as                  = it->second.as;
-                c[i].gfw                 = it->gfw;
+                c[i].gfw                 = it->second.gfw;
                 //c[i].skip                = 0;
                 c[i].phase               = NULL;
                 i++;
@@ -98,6 +99,47 @@ struct conc *cxxISolutionComp::cxxISolutionComp2conc(const std::set <cxxISolutio
         return(c);
 }
 
+void cxxISolutionComp::set_gfw() 
+{
+// return gfw
+  if (this->gfw > 0.0) return;
+// calculate gfw from as or from master species gfw
+  if (this->as != NULL)
+  {
+    /* use given chemical formula to calculate gfw */
+    double gfw;
+    if (compute_gfw (this->as, &gfw) == ERROR)
+    {
+      std::ostringstream oss;
+      oss << "Could not compute gfw, " <<  this->as;
+      error_msg(oss.str().c_str(), CONTINUE);
+      input_error++;
+      return;
+    }
+    if (this->description == "Alkalinity" && this->as == "CaCO3") 
+    {
+      gfw /= 2.;
+    }
+    this->gfw = gfw;
+    return;
+  }
+  /* use gfw of master species */
+  std::string str(this->description);
+  struct master *master_ptr = master_bsearch (str.c_str());
+  if (master_ptr != NULL)
+  {
+    /* use gfw for element redox state */
+    this->gfw = master_ptr->gfw;
+    return;
+  }
+  std::ostringstream oss;
+  oss << "Could not find gfw, " <<  this->description;
+  error_msg(oss.str().c_str(), CONTINUE);
+  input_error++;
+  return;
+}
+
+
 #ifdef SKIP
 cxxISolutionComp::STATUS_TYPE cxxISolutionComp::read(CParser& parser, cxxISolution& solution)
 {

ISolutionComp.h

@@ -3,6 +3,10 @@
 
 //#include "Parser.h"
 #include "Utils.h"
+#define EXTERNAL extern
+#include "global.h"
+#include "phqalloc.h"
+#include "phrqproto.h"
 
 #include <string>
 #include <map>     // std::map
@@ -26,24 +30,38 @@ public:
 
         void dump_xml(std::ostream& os, unsigned int indent = 0)const;
 
-        double get_input_conc()const {return this->input_conc;}
-        void set_input_conc(double input_conc) {this->input_conc = input_conc;}
-
-        std::string get_equation_name()const {return this->equation_name;}
-        void set_equation_name(char * equation_name) {this->equation_name = equation_name;}
-
         std::string get_description()const {return this->description;}
         void set_description(char * description) {this->description = description;}
 
+        double get_moles()const {return this->moles;}
+        void set_moles(double moles) {this->moles = moles;}
+
+        double get_input_conc()const {return this->input_conc;}
+        void set_input_conc(double input_conc) {this->input_conc = input_conc;}
+
         std::string get_units()const {return this->units;}
         void set_units(char * units) {this->units = units;}
 
+        std::string get_equation_name()const {return this->equation_name;}
+        void set_equation_name(char * equation_name) {this->equation_name = equation_name;}
+
+        double get_phase_si()const {return this->phase_si;}
+        void set_phase_si(int phase_si) {this->phase_si = phase_si;}
+
         int get_n_pe()const {return this->n_pe;}
         void set_n_pe(int n_pe) {this->n_pe = n_pe;}
 
+        char *get_as()const {return this->as;}
+        void set_as(char *as) {this->as = as;}
+
+        //double get_gfw()const {return this->gfw;}
+	double get_gfw()const {return this->gfw;};
+        void set_gfw(double gfw) {this->gfw = gfw;}
+	void set_gfw();
+
         bool operator<(const cxxISolutionComp& conc)const    { return ::strcmp(this->description, conc.description) < 0; }
 
-        static struct conc * cxxISolutionComp2conc(const std::set<cxxISolutionComp> &t );
+        static struct conc * cxxISolutionComp2conc(const std::map <char *, cxxISolutionComp, CHARSTAR_LESS> &t );
 
 private:
         char * description;

NameDouble.cxx

@@ -290,6 +290,18 @@ void cxxNameDouble::add(const cxxNameDouble &old, double factor)
         }
 
 }
+void cxxNameDouble::add(char * key, double total)
+//
+// add to total for a specified element
+//
+{
+  cxxNameDouble::iterator current = (*this).find(key);
+  if (current != (*this).end()) {
+    (*this)[key] = current->second + total;
+  } else {
+    (*this)[key] = total;
+  }
+}
 #ifdef USE_MPI
 #include "Dictionary.h"
 void cxxNameDouble::mpi_pack(std::vector<int>& ints, std::vector<double>& doubles) {

NameDouble.h

@@ -45,6 +45,7 @@ public:
         CParser::STATUS_TYPE read_raw(CParser& parser, std::istream::pos_type& pos);
 
         void add(const cxxNameDouble &old, double factor);
+	void add(char * key, double total);
 
         void insert(char *str, double d) {
 		(*this)[str] = d;

Parser.cxx

@@ -7,6 +7,7 @@
 
 #include "Parser.h"
 #include "Utils.h"
+#include "output.h"
 #include <algorithm>    // std::transform
 #include <map>          // std::map
 #include <cassert>      // assert
@@ -24,6 +25,8 @@ CParser::CParser(std::istream& input)
 {
         m_line_save.reserve(80);
         m_line.reserve(80);
+	echo_file = EO_ALL;
+	echo_stream = EO_NONE;
 }
 
 CParser::CParser(std::istream& input, std::ostream& output)
@@ -32,6 +35,8 @@ CParser::CParser(std::istream& input, std::ostream& output)
 {
         m_line_save.reserve(80);
         m_line.reserve(80);
+	echo_file = EO_ALL;
+	echo_stream = EO_NONE;
 }
 
 CParser::CParser(std::istream& input, std::ostream& output, std::ostream& error)
@@ -40,6 +45,8 @@ CParser::CParser(std::istream& input, std::ostream& output, std::ostream& error)
 {
         m_line_save.reserve(80);
         m_line.reserve(80);
+	echo_file = EO_ALL;
+	echo_stream = EO_NONE;
 }
 
 CParser::~CParser()
@@ -58,15 +65,50 @@ CParser::LINE_TYPE CParser::check_line(const std::string& str, bool allow_empty,
                 m_line_iss.seekg(0, std::ios_base::beg);
                 m_line_iss.clear();
 
-
+		// output for stream
-                if (true) // pr.echo_input == TRUE
+                switch (this->echo_stream) 
-                {
+		{
-                        if ((print && i != LT_EOF) || i == LT_KEYWORD)
+		case EO_NONE:
-                        {
+		  break;
-                                get_output() << "\t" << m_line_save << "\n";
+		case EO_ALL:
-                        }
+		  if (i != LT_EOF) 
-                }
+		  {
-
+		    std::ostringstream msg;
+		    msg << "\t" << m_line_save << "\n";
+		    get_output() << msg;
+		  }
+		  break;
+		case EO_KEWORDS:
+		  if (i == LT_KEYWORD) 
+		  {
+		    std::ostringstream msg;
+		    msg << "\t" << m_line_save << "\n";
+		    get_output() << msg;
+		  }
+		  break;
+		}
+		// output for file
+                switch (this->echo_file) 
+		{
+		case EO_NONE:
+		  break;
+		case EO_ALL:
+		  if (i != LT_EOF) 
+		  {
+		    std::ostringstream msg;
+		    msg << "\t" << m_line_save << "\n";
+		    output_msg(OUTPUT_MESSAGE, "%s", msg.str().c_str());
+		  }
+		  break;
+		case EO_KEWORDS:
+		  if (i == LT_KEYWORD) 
+		  {
+		    std::ostringstream msg;
+		    msg << "\t" << m_line_save << "\n";
+		    output_msg(OUTPUT_MESSAGE, "%s", msg.str().c_str());
+		  }
+		  break;
+		}
         } while (i == LT_EMPTY && allow_empty == false);
 
         // Check eof
@@ -507,7 +549,7 @@ int CParser::get_option(const std::vector<std::string>& opt_list, std::string::i
         //
         // Read line
         //
-        LINE_TYPE lt = check_line("get_option", false, true, true, false);
+        LINE_TYPE lt = check_line("get_option", false, true, true, true);
         if (lt == LT_EOF)
         {
                 j = OPT_EOF;
@@ -598,7 +640,7 @@ int CParser::get_option(const std::vector<std::string>& opt_list, std::istream::
         //
         // Read line
         //
-        LINE_TYPE lt = check_line("get_option", false, true, true, false);
+        LINE_TYPE lt = check_line("get_option", false, true, true, true);
         if (lt == LT_EOF)
         {
                 j = OPT_EOF;
@@ -678,6 +720,7 @@ int CParser::get_option(const std::vector<std::string>& opt_list, std::istream::
 
 int CParser::error_msg(const char *err_str, ONERROR_TYPE ot)
 {
+        ::error_msg(err_str, (int)ot);
         m_error_stream  << "ERROR: " << err_str << "\n";
         m_error_stream.flush();
 

Parser.h

@@ -1,6 +1,7 @@
 #if !defined(PARSER_H_INCLUDED)
 #define PARSER_H_INCLUDED
 
+extern int input_error;
 #include <string>   // std::string
 #include <map>      // std::map
 #include <vector>   // std::vector
@@ -64,12 +65,19 @@ public:
                 OT_CONTINUE     = 0,
                 OT_STOP         = 1
         };
+	
+	enum ECHO_OPTION {
+                EO_NONE         = 0,
+                EO_ALL          = 1,
+		EO_KEWORDS      = 2
+        };
 
-        enum STATUS_TYPE {
+	enum STATUS_TYPE {
                 PARSER_ERROR    = 0,
                 PARSER_OK       = 1
         };
 
+
         /**
                 Function gets a new line and checks for empty, eof, and keywords.
 
@@ -123,7 +131,7 @@ public:
 
         std::string& line()            {return m_line;}
         std::istringstream& get_iss()  {return m_line_iss;}
-        int incr_input_error()         {return ++m_input_error;}
+	int incr_input_error()         {++::input_error; return ++m_input_error;}
         std::ostream& get_output()     {return m_output_stream;}
         int get_input_error()          {return m_input_error;}
 
@@ -186,7 +194,12 @@ public:
         int error_msg(const std::ostringstream& err_str, ONERROR_TYPE stop)    {return error_msg(err_str.str().c_str(), stop);}
         int error_msg(const char *err_str, ONERROR_TYPE stop);
         int warning_msg(const char *err_str);
+        
+	void set_echo_file(ECHO_OPTION opt) {echo_file = opt;}
+	ECHO_OPTION get_echo_file() {return this->echo_file;};
 
+	void set_echo_stream(ECHO_OPTION opt) {echo_stream = opt;}
+	ECHO_OPTION get_echo_stream() {return this->echo_stream;};
 
         STATUS_TYPE parse_couple(std::string& token);
 
@@ -206,6 +219,8 @@ private:
         std::string          m_line_save;
         std::istringstream   m_line_iss;
         LINE_TYPE            m_line_type;
+	ECHO_OPTION          echo_stream;
+	ECHO_OPTION          echo_file;
 };
 
 #endif // PARSER_H_INCLUDED

ReadClass.cxx

@@ -17,6 +17,8 @@
 #include "phqalloc.h"
 #include "output.h"
 #include "phrqproto.h"
+#include <iostream>
+
 
 extern int check_line(const char *string, int allow_empty, int allow_eof, int allow_keyword,
 		      int print);
@@ -50,7 +52,7 @@ int read_solution_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("solution_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -119,7 +121,7 @@ int read_exchange_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("exchange_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("exchange_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -187,7 +189,7 @@ int read_surface_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("surface_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("surface_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -255,7 +257,7 @@ int read_equilibrium_phases_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("equilibrium_phases_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("equilibrium_phases_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -323,7 +325,7 @@ int read_kinetics_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("kinetics_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("kinetics_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -391,7 +393,7 @@ int read_solid_solutions_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -459,7 +461,7 @@ int read_gas_phase_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("gas_phase_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -527,7 +529,7 @@ int read_reaction_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("reaction_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -594,7 +596,7 @@ int read_mix_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("mix_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);
@@ -661,7 +663,7 @@ int read_temperature_raw (void)
  *   Read additonal lines
  */
 	for (;;) {
-		return_value = check_line("solid_solution_raw",TRUE,TRUE,TRUE,TRUE);
+		return_value = check_line("temperature_raw",TRUE,TRUE,TRUE,FALSE);
                                                /* empty, eof, keyword, print */
 		if (return_value == EOF || return_value == KEYWORD ) break;
 		keywordLines.append(line);

Solution.cxx

@@ -464,7 +464,7 @@ void cxxSolution::read_raw(CParser& parser)
                         opt = CParser::OPT_EOF;
                         parser.error_msg("Unknown input in SOLUTION_RAW keyword.", CParser::OT_CONTINUE);
                         parser.error_msg(parser.line().c_str(), CParser::OT_CONTINUE);
-                        break;
+                        continue;
 
                 case 0: // totals
                         if ( this->totals.read_raw(parser, next_char) != CParser::PARSER_OK) {

Utils.h

@@ -29,7 +29,7 @@ namespace Utilities {
 
         void squeeze_white(std::string& s_l);
 
-        void error_msg(const std::string&, const int stopflag);
+        //void error_msg(const std::string&, const int stopflag);
 
 	// operations on maps of entities (Solution, Exchange, ...)
         template<typename T>