Can read a solution result from Orchestra.cxx

Now will write routines to put results into unknown structures so that results can be printed.

git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@2166 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2007-07-30 15:07:20 +00:00
parent 1848045bd9
commit ea61d055f0
4 changed files with 224 additions and 29 deletions

View File

@ -30,6 +30,7 @@ units("mMol/kgw")
{
density = 1.0;
default_pe = -1;
pes = NULL;
}
cxxISolution::cxxISolution(struct solution *solution_ptr)
@ -419,7 +420,7 @@ void cxxISolution::dump_xml(std::ostream& os, unsigned int indent)const
}
#endif
void cxxISolution::ORCH_write_chemistry(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry(std::ostream &chemistry_dat)
{
this->ORCH_write_chemistry_water(chemistry_dat);
this->ORCH_write_chemistry_primary(chemistry_dat);
@ -430,7 +431,7 @@ void cxxISolution::ORCH_write_chemistry(std::ostringstream &chemistry_dat)
}
void cxxISolution::ORCH_write_chemistry_water(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_water(std::ostream &chemistry_dat)
{
//
// Write water entities
@ -488,7 +489,7 @@ void cxxISolution::ORCH_write_chemistry_water(std::ostringstream &chemistry_dat)
chemistry_dat << std::endl;
}
void cxxISolution::ORCH_write_chemistry_primary(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_primary(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* The primary species" << std::endl;
//
@ -535,7 +536,7 @@ void cxxISolution::ORCH_write_chemistry_primary(std::ostringstream &chemistry_da
chemistry_dat << std::endl;
}
void cxxISolution::ORCH_write_chemistry_total_O_H(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Define total hydrogen and oxygen" << std::endl;
// Define total hydrogen, total oxygen, and difference
@ -587,7 +588,7 @@ void cxxISolution::ORCH_write_chemistry_total_O_H(std::ostringstream &chemistry_
chemistry_dat << std::endl;
}
void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Alkalinity definitions" << std::endl;
//
@ -655,7 +656,7 @@ void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostringstream &chemistry
}
}
}
void cxxISolution::ORCH_write_chemistry_species(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_species(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Aqueous species" << std::endl;
//
@ -690,7 +691,7 @@ void cxxISolution::ORCH_write_chemistry_species(std::ostringstream &chemistry_da
}
chemistry_dat << std::endl;
}
void cxxISolution::ORCH_write_chemistry_minerals(std::ostringstream &chemistry_dat)
void cxxISolution::ORCH_write_chemistry_minerals(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Adjustments to mineral equilibrium" << std::endl;
//
@ -731,7 +732,7 @@ void cxxISolution::ORCH_write_chemistry_minerals(std::ostringstream &chemistry_d
}
}
}
void cxxISolution::ORCH_write_input(std::ostringstream &input_dat)
void cxxISolution::ORCH_write_input(std::ostream &input_dat)
{
@ -749,6 +750,8 @@ void cxxISolution::ORCH_write_input(std::ostringstream &input_dat)
//s_oss << "SOLUTION_RAW " << this->n_user << " " << this->description << std::endl;
//s_oss << "-temp " << this->tc << std::endl;
headings << "tempc\t";
data << this->tc << "\t";
//s_oss << "-pH " << this->ph << std::endl;
headings << "pH\t";
@ -847,18 +850,39 @@ void cxxISolution::ORCH_write_input(std::ostringstream &input_dat)
return;
}
void cxxISolution::ORCH_write_output(std::ostringstream &outstream)
void cxxISolution::ORCH_write_output(std::ostream &outstream)
{
outstream << "Var:";
outstream << "\tnr_iter";
//
// Write total concentrations in solution
// Serialize solution
//
outstream << "\tstart_solution";
//tc
outstream << "\ttempc";
//pH
outstream << "\tpH";
//pe
outstream << "\tpe";
//mu
outstream << "\tI";
//ah2o
outstream << "\tH2O.act";
//total_h;
outstream << "\ttotal_hydrogen";
//total_o;
outstream << "\ttotal_oxygen";
//cb
outstream << "\tchargebalance";
//mass_water;
outstream << "\tH2O.diss";
outstream << "\te-.diss";
//total_alkalinity;
outstream << "\tAlkalinity";
//orchestra master variables
outstream << "\tH+.diss";
std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::iterator iter = this->comps.begin();
for(; iter != this->comps.end(); ++iter)
outstream << "\te-.diss";
//totals
for(std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::iterator iter = this->comps.begin(); iter != this->comps.end(); ++iter)
{
std::string name(iter->second.get_description());
if (name == "H(1)" || name == "E" || name == "Alkalinity") continue;
@ -871,8 +895,21 @@ void cxxISolution::ORCH_write_output(std::ostringstream &outstream)
assert(s_ptr != NULL);
outstream << "\t" << s_ptr->name << ".diss";
}
outstream << "\tchargebalance";
outstream<< "\tI";
outstream << "\tend_totals";
for(std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::iterator iter = this->comps.begin(); iter != this->comps.end(); ++iter)
{
std::string name(iter->second.get_description());
if (name == "H(1)" || name == "E" || name == "Alkalinity") continue;
struct element *elt;
char *element_name = string_hsave(name.c_str());
elt = element_store(element_name);
assert(elt != NULL);
struct species *s_ptr;
s_ptr = elt->master->s;
assert(s_ptr != NULL);
outstream << "\t" << s_ptr->name << ".act";
}
outstream << "\tend_master_activities";
//
// Write all species activities and concentrations
//
@ -883,6 +920,89 @@ void cxxISolution::ORCH_write_output(std::ostringstream &outstream)
std::replace(name.begin(), name.end(), ')', ']');
outstream << "\t" << name.c_str() << ".act" << "\t" << name.c_str() << ".con";
}
outstream << "\tend_species";
outstream << std::endl;
return;
}
/* ---------------------------------------------------------------------- */
int
initial_solutions (int print)
/* ---------------------------------------------------------------------- */
{
/*
* Go through list of solutions, make initial solution calculations
* for any marked "new".
*/
int i, converge, converge1;
int n, last, n_user, print1;
char token[2 * MAX_LENGTH];
state = INITIAL_SOLUTION;
set_use ();
for (n = 0; n < count_solution; n++)
{
initial_solution_isotopes = FALSE;
if (solution[n] != NULL && solution[n]->new_def == TRUE)
{
if (print1 == TRUE && print == TRUE)
{
dup_print ("Beginning of initial solution calculations.", TRUE);
print1 = FALSE;
}
if (print == TRUE)
{
sprintf (token, "Initial solution %d.\t%.350s", solution[n]->n_user,
solution[n]->description);
dup_print (token, FALSE);
}
use.solution_ptr = solution[n];
prep ();
k_temp (solution[n]->tc);
set (TRUE);
converge = model ();
if (converge == ERROR && diagonal_scale == FALSE)
{
diagonal_scale = TRUE;
set (TRUE);
converge = model ();
diagonal_scale = FALSE;
}
converge1 = check_residuals ();
sum_species ();
add_isotopes (solution[n]);
punch_all ();
print_all ();
/* free_model_allocs(); */
if (converge == ERROR || converge1 == ERROR)
{
error_msg ("Model failed to converge for initial solution.", STOP);
}
n_user = solution[n]->n_user;
last = solution[n]->n_user_end;
/* copy isotope data */
if (solution[n]->count_isotopes > 0)
{
count_isotopes_x = solution[n]->count_isotopes;
isotopes_x =
(struct isotope *) PHRQ_realloc (isotopes_x,
(size_t) count_isotopes_x *
sizeof (struct isotope));
if (isotopes_x == NULL)
malloc_error ();
memcpy (isotopes_x, solution[n]->isotopes,
(size_t) count_isotopes_x * sizeof (struct isotope));
}
else
{
count_isotopes_x = 0;
}
xsolution_save (n_user);
for (i = n_user + 1; i <= last; i++)
{
solution_duplicate (n_user, i);
}
}
}
initial_solution_isotopes = FALSE;
return (OK);
}

View File

@ -41,9 +41,9 @@ public:
//void dump_xml(std::ostream& os, unsigned int indent = 0)const;
void ConvertUnits();
void ORCH_write_chemistry(std::ostringstream &chemistry_dat);
void ORCH_write_input(std::ostringstream &input_dat);
void ORCH_write_output(std::ostringstream &input_dat);
void ORCH_write_chemistry(std::ostream &chemistry_dat);
void ORCH_write_input(std::ostream &input_dat);
void ORCH_write_output(std::ostream &input_dat);
void print();
protected:
friend class cxxISolutionComp; // for this->pe access
@ -57,12 +57,12 @@ public:
//static std::map<int, cxxISolution>& map;
private:
void ORCH_write_chemistry_water(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_primary(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_total_O_H(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_alkalinity(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_species(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_minerals(std::ostringstream &chemistry_dat);
void ORCH_write_chemistry_water(std::ostream &chemistry_dat);
void ORCH_write_chemistry_primary(std::ostream &chemistry_dat);
void ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat);
void ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat);
void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
void ORCH_write_chemistry_minerals(std::ostream &chemistry_dat);
};
#endif // !defined(ISOLUTION_H_INCLUDED)

View File

@ -343,8 +343,8 @@ void cxxSolution::dump_xml(std::ostream& s_oss, unsigned int indent)const
s_oss << indent1;
s_oss << "soln_total_o=\"" << this->total_o << "\"" << std::endl;
s_oss << indent1;
s_oss << indent1;
s_oss << "soln_cb=\"" << this->cb << "\"" << std::endl;
s_oss << indent1;
@ -446,7 +446,7 @@ void cxxSolution::dump_raw(std::ostream& s_oss, unsigned int indent)const
s_oss << indent1;
s_oss << "-total_o " << this->total_o << std::endl;
// new identifier
// new identifier
s_oss << indent1;
s_oss << "-cb " << this->cb << std::endl;
@ -819,7 +819,8 @@ void cxxSolution::read_raw(CParser& parser)
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break;
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break;
}
// all members must be defined
if (tc_defined == false) {
@ -862,7 +863,7 @@ void cxxSolution::read_raw(CParser& parser)
parser.incr_input_error();
parser.error_msg("Total alkalinity not defined for SOLUTION_RAW input.", CParser::OT_CONTINUE);
}
return;
return;
}
void cxxSolution::zero()
@ -963,6 +964,79 @@ double cxxSolution::get_master_activity(char *string)const
return(it->second);
}
}
void cxxSolution::read_orchestra(std::vector <std::pair <std::string, double> > output_vector, std::vector < std::pair < std::string, double > >::iterator &it)
{
this->tc = it->second; it++;
this->ph = it->second; it++;
this->pe = it->second; it++;
this->mu = it->second; it++;
this->ah2o = it->second; it++;
this->total_h = it->second; it++;
this->total_o = it->second; it++;
this->cb = it->second; it++;
this->mass_water = it->second * gfw_water; it++;
this->total_alkalinity = it->second; it++;
it++; //orch total H+
it++; //orch total e-
//cxxNameDouble totals;
char token[MAX_LENGTH];
while (it->first.compare("end_totals") != 0)
{
strcpy(token, it->first.c_str());
replace(".diss","",token);
struct species *s_ptr = s_search(token);
if (s_ptr == NULL) error_msg("Species not found in orchestra read", STOP);
if (s_ptr->secondary != NULL)
{
this->totals[s_ptr->secondary->elt->name] = it->second;
} else if (s_ptr->primary != NULL)
{
this->totals[s_ptr->primary->elt->name] = it->second;
} else
{
error_msg("Species not secondary or primary master species in orchestra read", STOP);
}
it++;
}
//cxxNameDouble master_activity;
it++;
while (it->first.compare("end_master_activities") != 0)
{
strcpy(token, it->first.c_str());
replace(".act","",token);
struct species *s_ptr = s_search(token);
if (s_ptr == NULL) error_msg("Species not found in orchestra read", STOP);
if (s_ptr->secondary != NULL)
{
this->master_activity[s_ptr->secondary->elt->name] = log10(it->second);
} else if (s_ptr->primary != NULL)
{
this->master_activity[s_ptr->primary->elt->name] = log10(it->second);
} else
{
error_msg("Species not secondary or primary master species in orchestra read", STOP);
}
it++;
}
//cxxNameDouble species_gamma;
//cxxSolutionIsotopeList isotopes;
//
// Also process aqueous species data
//
it++;
while (it->first.compare("end_species") != 0)
{
strcpy(token, it->first.c_str());
replace(".act","",token);
while (replace("[","(", token));
while (replace("]",")", token));
struct species *s_ptr = s_search(token);
if (s_ptr == NULL) error_msg("Species not found in orchestra read", STOP);
s_ptr->la = log10(it->second); it++;
s_ptr->lm = log10(it->second); it++;
s_ptr->lg = s_ptr->la - s_ptr->lm;
}
}
#ifdef USE_MPI
#include <mpi.h>
/* ---------------------------------------------------------------------- */

View File

@ -84,6 +84,7 @@ public:
void read_raw(CParser& parser);
void multiply(double extensive);
void read_orchestra(std::vector <std::pair <std::string, double>> output_vector, std::vector < std::pair < std::string, double > >::iterator &it);
#ifdef USE_MPI