max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 16:44:49 +01:00
Code Issues Packages Projects Releases Wiki Activity

Removed ORCHESTRA

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/branches/ErrorHandling@5658 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2011-09-28 22:30:46 +00:00
parent 4d728afcbe
commit 6acf05a1a0
13 changed files with 9 additions and 1101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
GasPhase.h
View File

@ -9,6 +9,7 @@
#include "NumKeyword.h"
#include "NameDouble.h"
#include "phreeqc_class.h"
class cxxMix;
class cxxGasPhase:public cxxNumKeyword

591
ISolution.cxx
View File

@ -20,9 +20,6 @@
#include "phrqproto.h"
#include "output.h"
#ifdef ORCHESTRA
extern void ORCH_write_chemistry_species(std::ostream & chemistry_dat);
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
@ -460,591 +457,3 @@ cxxISolution::dump_xml(std::ostream & os, unsigned int indent) const const
os << "</solution>" << "\n";
}
#endif
#ifdef ORCHESTRA
void
cxxISolution::ORCH_write_chemistry(std::ostream & chemistry_dat)
{
this->ORCH_write_chemistry_water(chemistry_dat);
this->ORCH_write_chemistry_primary(chemistry_dat);
this->ORCH_write_chemistry_total_O_H(chemistry_dat);
this->ORCH_write_chemistry_alkalinity(chemistry_dat);
ORCH_write_chemistry_species(chemistry_dat);
this->ORCH_write_chemistry_minerals(chemistry_dat);
}
void
cxxISolution::ORCH_write_chemistry_water(std::ostream & chemistry_dat)
{
//
// Write water entities
//
chemistry_dat << std::endl << "//********* The water entities" << std::
endl;
// e-
chemistry_dat << "@entity(e-, diss, 0)" << std::endl;
if (this->comps.find("E") == this->comps.end()
|| this->comps.find("E")->second.get_equation_name() == NULL)
{
// fixed pe
chemistry_dat << "@Calc: (1, \"e-.act = 10^-pe\")" << std::endl;
}
else if (this->comps.find("E")->second.get_equation_name() ==
string_hsave("charge"))
{
// charge balance
chemistry_dat << "@Calc: (1, \"pe = -log10({e-.act})\")" << std::endl;
chemistry_dat <<
"@solve (e-.act, 1e-6, log, 1, chargebalance, 1e-14)" << std::
endl;
}
else
{
// adjust to equilbirium with a phase
chemistry_dat << "@Calc: (1, \"pe = -log10({e-.act})\")" << std::endl;
int n;
struct phase *phase_ptr =
phase_bsearch(this->comps.find("E)")->second.get_equation_name(),
&n, FALSE);
assert(phase_ptr != NULL);
std::string phase_name(phase_ptr->name);
std::replace(phase_name.begin(), phase_name.end(), '(', '[');
std::replace(phase_name.begin(), phase_name.end(), ')', ']');
chemistry_dat << "@solve (e-.act, 1e-6, log, 1, " << phase_name.
c_str() << ".si_raw, 1e-9)" << std::endl;
}
// H+
if (this->comps.find("H(1)") == this->comps.end()
|| this->comps.find("H(1)")->second.get_equation_name() == NULL)
{
// fixed pH
chemistry_dat << "@Calc: (1, \"H+.act = 10^-pH\")" << std::endl;
}
else if (this->comps.find("H(1)")->second.get_equation_name() ==
string_hsave("charge"))
{
// charge balance
chemistry_dat << "@Calc: (1, \"pH = -log10({H+.act})\")" << std::endl;
chemistry_dat <<
"@solve (H+.act, 1e-6, log, 1, chargebalance, 1e-14)" << std::
endl;
}
else
{
// adjust to equilbirium with a phase
chemistry_dat << "@Calc: (1, \"pH = -log10({H+.act})\")" << std::endl;
int n;
struct phase *phase_ptr =
phase_bsearch(this->comps.find("H(1)")->second.
get_equation_name(), &n, FALSE);
assert(phase_ptr != NULL);
chemistry_dat << "@solve (H+.act, 1e-6, log, 1, " << phase_ptr->
name << ".si_raw, 1e-9)" << std::endl;
}
// H2O
chemistry_dat << "@entity(" << s_h2o->
name << ", diss, 55.506)" << std::endl;
chemistry_dat << std::endl;
}
void
cxxISolution::ORCH_write_chemistry_primary(std::ostream & chemistry_dat)
{
chemistry_dat << std::endl << "//********* The primary species" << std::
endl;
//
// Write other master species definitions, i.e. primary entities
//
std::map < char *, cxxISolutionComp, CHARSTAR_LESS >::iterator iter =
this->comps.begin();
chemistry_dat << "@species(H+, 1)" << std::endl;
for (; 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);
chemistry_dat << "@species(" << s_ptr->name << ", " << s_ptr->
z << ")" << std::endl;
if (iter->second.get_equation_name() == NULL)
{
// regular mole balance
chemistry_dat << "@primary_entity(" << s_ptr->
name << ", 1e-9, diss, 1.0e-9)" << std::endl;
}
else
{
std::string eqn(iter->second.get_equation_name());
if (eqn == "charge")
{
// charge balance
chemistry_dat << "@solve (" << s_ptr->
name << ".act, 1e-6, log, 1, chargebalance, 1e-9)" <<
std::endl;
}
else
{
// adjust to phase equilibrium
int n;
struct phase *phase_ptr =
phase_bsearch(eqn.c_str(), &n, FALSE);
assert(phase_ptr != NULL);
std::string phase_name(phase_ptr->name);
std::replace(phase_name.begin(), phase_name.end(), '(', '[');
std::replace(phase_name.begin(), phase_name.end(), ')', ']');
chemistry_dat << "@solve (" << s_ptr->
name << ".act, 1e-6, log, 1, " << phase_name <<
".si_raw, 1e-9)" << std::endl;
}
}
}
chemistry_dat << std::endl;
}
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
//chemistry_dat << "@var: total_diff 0" << std::endl;
//chemistry_dat << "@calc: (5, \"total_diff = total_hydrogen - 2*total_oxygen" << "\")" << std::endl;
// Write H total equation
chemistry_dat << "@var: total_hydrogen 0" << std::endl;
chemistry_dat << "@calc: (5, \"total_hydrogen = 0";
int i;
for (i = 0; i < count_s_x; i++)
{
// write in terms of orchestra components
if (s_x[i]->primary != NULL
|| (s_x[i]->secondary != NULL && s_x[i]->secondary->in == TRUE))
{
if (s_x[i]->h != 0)
{
chemistry_dat << "+";
if (s_x[i]->h != 1)
{
chemistry_dat << s_x[i]->h << "*";
}
chemistry_dat << "{" << s_x[i]->name << ".diss}";
}
}
}
chemistry_dat << "\")" << std::endl;
// Write O total equation
chemistry_dat << "@var: total_oxygen 0" << std::endl;
chemistry_dat << "@calc: (5, \"total_oxygen = 0";
for (i = 0; i < count_s_x; i++)
{
if (s_x[i]->o != 0)
{
// write in terms of orchestra components
if (s_x[i]->primary != NULL
|| (s_x[i]->secondary != NULL
&& s_x[i]->secondary->in == TRUE))
{
chemistry_dat << "+";
if (s_x[i]->o != 1)
{
chemistry_dat << s_x[i]->o << "*";
}
chemistry_dat << "{" << s_x[i]->name << ".diss}";
}
}
}
chemistry_dat << "\")" << std::endl;
chemistry_dat << std::endl;
}
void
cxxISolution::ORCH_write_chemistry_alkalinity(std::ostream & chemistry_dat)
{
chemistry_dat << std::
endl << "//********* Alkalinity definitions" << std::endl;
//
// Define alkalinity
//
chemistry_dat << "@Var: Alkalinity 0" << std::endl;
chemistry_dat << "@calc: (5, \"Alkalinity = 0";
for (int i = 0; i < count_s_x; i++)
{
if (s_x[i]->alk == 0)
continue;
std::string name(s_x[i]->name);
std::replace(name.begin(), name.end(), '(', '[');
std::replace(name.begin(), name.end(), ')', ']');
if (s_x[i]->alk < 0)
{
if (s_x[i]->alk != -1)
{
chemistry_dat << s_x[i]->alk << "*{" << name << ".con}";
}
else
{
chemistry_dat << "-{" << name << ".con}";
}
}
else if (s_x[i]->alk > 0)
{
if (s_x[i]->alk != 1)
{
chemistry_dat << "+" << s_x[i]->
alk << "*{" << name << ".con}";
}
else
{
chemistry_dat << "+{" << name << ".con}";
}
}
}
chemistry_dat << "\")" << std::endl;
//
// Alkalinity (or pH) equation
//
std::map < char *, cxxISolutionComp, CHARSTAR_LESS >::iterator iter =
this->comps.begin();
if ((iter = this->comps.find("Alkalinity")) != this->comps.end())
{
if ((this->comps.find("C(4)") != this->comps.end())
|| (this->comps.find("C") != this->comps.end()))
{
if (this->comps.find("H(1)") != this->comps.end())
{
std::ostringstream oss;
oss <<
"pH can not be adjusted to charge balance or phase equilibrium when specifying C or C(4) and Alkalinty.";
error_msg(oss.str().c_str(), CONTINUE);
//P_INSTANCE_POINTER input_error++;
}
chemistry_dat << "@solve (pH, 1e-6, lin, 1, Alkalinity, 7)" <<
std::endl;
chemistry_dat << std::endl;
}
else
{
struct master *master_ptr = master_bsearch("Alkalinity");
if (master_ptr == NULL)
{
std::ostringstream oss;
oss << "Could not find Alkalinity definition in database.";
error_msg(oss.str().c_str(), CONTINUE);
input_error++;
}
chemistry_dat << "@species(" << master_ptr->s->
name << ", " << master_ptr->s->z << ")" << std::endl;
chemistry_dat << "@solve (" << master_ptr->s->
name << ".act, 1e-6, log, 1, Alkalinity, 1e-9)" << std::endl;
chemistry_dat << std::endl;
}
}
}
void
cxxISolution::ORCH_write_chemistry_minerals(std::ostream & chemistry_dat)
{
chemistry_dat << std::
endl << "//********* Adjustments to mineral equilibrium" << std::endl;
//
// Write minerals
//
std::map < char *, cxxISolutionComp, CHARSTAR_LESS >::iterator iter =
this->comps.begin();
for (iter = this->comps.begin(); iter != this->comps.end(); ++iter)
{
if (iter->second.get_equation_name() != NULL)
{
std::string name(iter->second.get_equation_name());
if (name != "charge")
{
struct phase *phase_ptr;
int n;
phase_ptr = phase_bsearch(name.c_str(), &n, FALSE);
assert(phase_ptr != NULL);
std::string phase_name(phase_ptr->name);
std::replace(phase_name.begin(), phase_name.end(), '(', '[');
std::replace(phase_name.begin(), phase_name.end(), ')', ']');
chemistry_dat << "@si_mineral(" << phase_name << ")" << std::
endl;
chemistry_dat << "@reaction(" << phase_name << ", " <<
pow(10.0, -phase_ptr->lk);
struct rxn_token *next_token = phase_ptr->rxn_x->token;
next_token++;
while (next_token->s != NULL || next_token->name != NULL)
{
chemistry_dat << ", " << next_token->coef;
if (next_token->s != NULL)
{
chemistry_dat << ", " << next_token->s->name;
}
else
{
chemistry_dat << ", " << next_token->name;
}
next_token++;
}
chemistry_dat << ")" << std::endl;
}
}
}
//chemistry_dat << "@mineral(Quartz)" << std::endl;
//chemistry_dat << "@sreaction(Quartz, 10139.1138573668, -2.0, H2O, 1.0, H4SiO4)" << std::endl;
}
void
cxxISolution::ORCH_write_input(std::ostream & input_dat)
{
//
// Write orchestra input file info
//
std::ostringstream headings, data;
data.precision(DBL_DIG - 1);
headings << "var: ";
data << "data: ";
// Solution element and attributes
//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";
data << this->ph << "\t";
//s_oss << "-pe " << this->pe << std::endl;
headings << "pe\t";
data << this->pe << "\t";
//s_oss << "-mu " << this->mu << std::endl;
//s_oss << "-ah2o " << this->ah2o << std::endl;
headings << "H2O.act\t";
data << 1 << "\t";
//s_oss << "-total_h " << this->total_h << std::endl;
//s_oss << "-total_o " << this->total_o << std::endl;
//s_oss << "-cb " << this->cb << std::endl;
//s_oss << "-mass_water " << this->mass_water << std::endl;
//s_oss << "-total_alkalinity " << this->total_alkalinity << std::endl;
// soln_total conc structures
//this->totals.dump_raw(s_oss, indent + 2);
//this->totals.write_orchestra(headings, s_oss);
std::map < char *, cxxISolutionComp, CHARSTAR_LESS >::iterator iter =
this->comps.begin();
for (; iter != this->comps.end(); ++iter)
{
std::string master_name;
struct master *master_ptr;
master_ptr = master_bsearch(iter->first);
assert(master_ptr != NULL);
std::string ename(iter->first);
double coef = master_ptr->coef;
if (master_ptr->coef == 0)
{
coef = 1;
}
if (ename != "Alkalinity")
{
ename = master_ptr->s->name;
ename.append(".diss");
}
if (iter->second.get_equation_name() == NULL)
{
headings << ename << "\t";
data << (this->totals.find(iter->first))->second / coef << "\t";
}
else
{
std::string name(iter->second.get_equation_name());
if (name == "charge")
{
headings << ename << "\t";
data << (this->totals.find(iter->first))->second /
coef << "\t";
}
else
{
int n;
struct phase *phase_ptr =
phase_bsearch(name.c_str(), &n, TRUE);
assert(phase_ptr != NULL);
std::string phase_name(phase_ptr->name);
std::replace(phase_name.begin(), phase_name.end(), '(', '[');
std::replace(phase_name.begin(), phase_name.end(), ')', ']');
headings << phase_name << ".si_raw" << "\t";
data << iter->second.get_phase_si() << "\t";
}
}
// activity estimate
if (ename == "Alkalinity")
{
if ((this->comps.find("C") == this->comps.end())
&& (this->comps.find("C(4)") == this->comps.end()))
{
headings << master_ptr->s->name << ".act\t";
data << iter->second.get_moles() * 1e-3 << "\t";
}
}
else
{
headings << master_ptr->s->name << ".act\t";
data << iter->second.get_moles() * 1e-3 << "\t";
}
}
// Isotopes
//s_oss << "-Isotopes" << std::endl;
/*
{
for (std::list<cxxSolutionIsotope>::const_iterator it = this->isotopes.begin(); it != isotopes.end(); ++it) {
it->dump_raw(s_oss, indent + 2);
}
}
*/
// Write data to string
input_dat << headings.str() << std::endl;
input_dat << data.str() << std::endl;
return;
}
void
cxxISolution::ORCH_write_output_vars(std::ostream & outstream)
{
outstream << "Var:";
outstream << "\tnr_iter";
//
// 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.con";
//total_alkalinity;
outstream << "\tAlkalinity";
//orchestra master variables
outstream << "\tH+.diss";
outstream << "\te-.diss";
outstream << "\tH2O.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;
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 << ".diss";
}
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
//
int i;
for (i = 0; i < count_s_x; i++)
{
std::string name(s_x[i]->name);
std::replace(name.begin(), name.end(), '(', '[');
std::replace(name.begin(), name.end(), ')', ']');
outstream << "\t" << name.c_str() << ".act" << "\t" << name.
c_str() << ".con";
}
outstream << "\tend_species";
outstream << std::endl;
return;
}
void
cxxISolution::ORCH_write(std::ostream & chemistry_dat,
std::ostream & input_dat, std::ostream & output_dat)
{
//
// Write orchestra chemistry file definition
//
chemistry_dat << std::
endl << "// Write ORCHESTRA chemistry definitions" << std::endl;
// mark for Orchestra include
chemistry_dat << std::endl << "@class: species_reactions () {" << std::
endl;
this->ORCH_write_chemistry(chemistry_dat);
// end orchestra include block
chemistry_dat << std::endl << "}" << std::endl;
//
// Write orchestra input file definition
//
input_dat << std::endl << "@class: input_file_data () {" << std::endl;
this->ORCH_write_input(input_dat);
input_dat << std::endl << "}" << std::endl;
//
// Write orchestra output file definition
//
output_dat << std::endl << "Output_every: 1" << std::endl;
this->ORCH_write_output_vars(output_dat);
//write data to stderr
//std::cerr << chemistry_dat.str() << input_dat.str() << output_dat.str();
}
#endif

1
NameDouble.h
View File

@ -7,7 +7,6 @@
#include <string> // std::string
#include <list> // std::list
#include <vector> // std::vector
#include "Phreeqc_class.h"
class Phreeqc;
#include "Parser.h"

92
PPassemblage.cxx
View File

@ -342,95 +342,3 @@ cxxPPassemblage::add(const cxxPPassemblage & addee, double extensive)
this->eltList.add_extensive(addee.eltList, extensive);
}
#ifdef ORCHESTRA
void
cxxPPassemblage::ORCH_write_chemistry(std::ostream & chemistry_dat)
{
chemistry_dat << std::endl << "//********* Mineral equilibria" << std::
endl;
//
// Write minerals
//
for (std::list < cxxPPassemblageComp >::iterator it =
this->ppAssemblageComps.begin(); it != this->ppAssemblageComps.end();
++it)
{
struct phase *phase_ptr;
int n;
phase_ptr = phase_bsearch(it->get_name(), &n, FALSE);
assert(phase_ptr != NULL);
std::string phase_name(phase_ptr->name);
std::replace(phase_name.begin(), phase_name.end(), '(', '[');
std::replace(phase_name.begin(), phase_name.end(), ')', ']');
chemistry_dat << "@mineral(" << phase_name << ")" << std::endl;
chemistry_dat << "@reaction(" << phase_name << ", " << pow(10.0,
-phase_ptr->
lk);
struct rxn_token *next_token = phase_ptr->rxn_x->token;
next_token++;
while (next_token->s != NULL || next_token->name != NULL)
{
chemistry_dat << ", " << next_token->coef;
if (next_token->s != NULL)
{
chemistry_dat << ", " << next_token->s->name;
}
else
{
chemistry_dat << ", " << next_token->name;
}
next_token++;
}
chemistry_dat << ")" << std::endl;
}
}
void
cxxPPassemblage::ORCH_write_output_vars(std::ostream & outstream)
{
//
// Write ppassemblage output variables
//
outstream << "\tstart_ppassemblage";
for (std::list < cxxPPassemblageComp >::iterator it =
this->ppAssemblageComps.begin(); it != this->ppAssemblageComps.end();
++it)
{
outstream << "\t" << it->get_name() << ".min";
}
outstream << "\tend_ppassemblage";
}
void
cxxPPassemblage::ORCH_read(std::vector < std::pair < std::string,
double >>output_vector,
std::vector < std::pair < std::string,
double >>::iterator & it)
{
while (it->first.compare("end_ppassemblage") != 0)
{
it++;
}
}
void
cxxPPassemblage::ORCH_store_global(std::map < std::string, double >output_map)
{
int i;
std::map < std::string, double >::iterator it;
for (i = 0; i < count_unknowns; i++)
{
// MB, ALK, CB, SOLUTION_PHASE_BOUNDARY, MU, AH2O
switch (x[i]->type)
{
case PP:
std::string name(x[i]->phase->name);
name.append(".min");
it = output_map.find(name);
assert(it != output_map.end());
x[i]->moles = it->second;
break;
}
}
}
#endif

9
PPassemblage.h
View File

@ -43,15 +43,6 @@ class cxxPPassemblage:public cxxNumKeyword
#endif
void totalize(PHREEQC_PTR_ARG);
#ifdef ORCHESTRA
void ORCH_write_chemistry(std::ostream & chemistry_dat);
void ORCH_write_output_vars(std::ostream & outstream);
void ORCH_read(std::vector < std::pair < std::string,
double >>output_vector,
std::vector < std::pair < std::string,
double >>::iterator & it);
void ORCH_store_global(std::map < std::string, double >output_map);
#endif
private:
void add(const cxxPPassemblage & addee, double extensive);

1
PPassemblageComp.h
View File

@ -8,6 +8,7 @@
#include <vector> // std::vector
#include "NameDouble.h"
#include "phreeqc_class.h"
class cxxPPassemblageComp: public PHRQ_base
{

2
SSassemblage.h
View File

@ -10,6 +10,8 @@
#include "NumKeyword.h"
#include "NameDouble.h"
#include "phreeqc_class.h"
class cxxSSassemblageSS;
//#include "cxxMix.h"

1
SSassemblageSS.h
View File

@ -7,6 +7,7 @@
#include <list> // std::list
#include <vector> // std::vector
#include "NameDouble.h"
#include "phreeqc_class.h"
class cxxSSassemblageSS: public PHRQ_base
{

9
Solution.h
View File

@ -11,6 +11,7 @@
#include "NameDouble.h"
#include "PHRQ_base.h"
#include "PHRQ_io.h"
#include "phreeqc_class.h"
class cxxMix;
class cxxSolution:public cxxNumKeyword
@ -151,14 +152,6 @@ class cxxSolution:public cxxNumKeyword
void read_raw(CParser & parser, bool check = true);
void multiply(double extensive);
void modify_activities(PHREEQC_PTR_ARG_COMMA const cxxSolution & original);
#ifdef ORCHESTRA
void ORCH_write(std::ostream & headings, std::ostream & input_data) const;
void ORCH_read(std::vector < std::pair < std::string,
double >>output_vector,
std::vector < std::pair < std::string,
double >>::iterator & it);
void ORCH_store_global(std::map < std::string, double >output_map);
#endif
#ifdef USE_MPI
void mpi_pack(std::vector < int >&ints, std::vector < double >&doubles);

16
StorageBin.cxx
View File

@ -922,7 +922,7 @@ cxxStorageBin::add(struct system *system_ptr)
void
cxxStorageBin::cxxStorageBin2phreeqc(PHREEQC_PTR_ARG_COMMA int n)
//
// copy data fromphreeqc storage to storage bin
// copy data from phreeqc storage to storage bin
//
{
//std::ostringstream oss;
@ -1047,7 +1047,7 @@ cxxStorageBin::cxxStorageBin2phreeqc(PHREEQC_PTR_ARG_COMMA int n)
void
cxxStorageBin::phreeqc2cxxStorageBin(PHREEQC_PTR_ARG_COMMA int n)
//
// copy data fromphreeqc storage to storage bin
// copy data from phreeqc storage to storage bin
//
{
int pos;
@ -1645,18 +1645,6 @@ cxxStorageBin::mix_cxxExchange(cxxMix & mixmap)
}
#endif
#ifdef ORCHESTRA
void
cxxStorageBin::ORCH_write(std::ostream & chemistry_dat,
std::ostream & input_dat, std::ostream & output_dat)
{
//std::ostringstream oss;
// Liter concentrations
this->system.ORCH_write(chemistry_dat, input_dat, output_dat);
}
#endif
cxxSystem &
cxxStorageBin::getSystem(void)
{

4
StorageBin.h
View File

@ -108,10 +108,6 @@ class cxxStorageBin: public PHRQ_base
#ifdef USE_MPI
void mpi_send(int n, int task_number);
void mpi_recv(int task_number);
#endif
#ifdef ORCHESTRA
void ORCH_write(std::ostream & chemistry_dat, std::ostream & input_dat,
std::ostream & output_dat);
#endif
protected:
// Tidied classes

370
System.cxx
View File

@ -20,8 +20,6 @@
#include "Reaction.h"
#include "Temperature.h"
extern void ORCH_write_chemistry_species(std::ostream & chemistry_dat);
cxxSystem::cxxSystem(PHRQ_io *io)
:
PHRQ_base(io)
@ -96,373 +94,5 @@ cxxSystem::totalize(PHREEQC_PTR_ARG)
this->ssassemblage->totalize(P_INSTANCE);
this->totals.add_extensive(this->surface->Get_totals(), 1.0);
}
//Need to handle the following 3 reactions:
//
//this->kinetics = NULL;
//this->mix = NULL;
//this->reaction = NULL;
//this->temperature = NULL;
//this->totals.dump_raw(std::cerr, 1);
return;
}
#ifdef ORCHESTRA
void
cxxSystem::ORCH_write(std::ostream & chemistry_dat, std::ostream & input_dat,
std::ostream & output_dat)
{
//
// chemistry file
//
this->totalize();
// calculate orchestra components from PHREEQC totals
this->ORCH_components();
// write definitions for chemistry
chemistry_dat << std::
endl << "@class: species_reactions () {" << std::endl;
this->ORCH_write_chemistry_water(chemistry_dat);
this->ORCH_write_chemistry_primary(chemistry_dat);
this->ORCH_write_chemistry_total_O_H(chemistry_dat);
ORCH_write_chemistry_species(chemistry_dat);
// add definitions for pure phases
if (this->ppassemblage != NULL)
this->ppassemblage->ORCH_write_chemistry(chemistry_dat);
// add definitions for other PHREEQC entities here
// finish up
chemistry_dat << std::endl << "}" << std::endl;
//
// input file
//
input_dat << std::endl << "@class: input_file_data () {" << std::endl;
this->ORCH_write_input(input_dat);
input_dat << std::endl << "}" << std::endl;
//
// output file
//
output_dat << std::endl << "Output_every: 1" << std::endl;
output_dat << "Var:";
this->ORCH_write_output_vars(output_dat);
// add definitions for pure phases
if (this->ppassemblage != NULL)
this->ppassemblage->ORCH_write_output_vars(output_dat);
//finish up
output_dat << std::endl;
};
void
cxxSystem::ORCH_write_chemistry_water(std::ostream & chemistry_dat)
{
//
// Write water entities
//
chemistry_dat << std::
endl << "//********* The water entities" << std::endl;
// e-, total hydrogen
chemistry_dat << "@entity(e-, diss, 0)" << std::endl;
chemistry_dat << "@Calc: (1, \"e-.act = 10.^(-pe)\")" << std::endl;
chemistry_dat << "@solve (pe, 1e-6, lin, 1, e-.liter, 1e-14)" << std::
endl;
// H+, charge balance
chemistry_dat << "@Calc: (1, \"H+.act = 10.^(-pH)\")" << std::endl;
chemistry_dat << "@solve (pH, 1e-6, lin, 1, H+.liter, 1e-14)" << std::
endl;
// H2O
chemistry_dat << "@entity(" << s_h2o->
name << ", diss, 55.506)" << std::endl;
chemistry_dat << std::endl;
}
void
cxxSystem::ORCH_write_chemistry_primary(std::ostream & chemistry_dat)
{
chemistry_dat << std::endl << "//********* The primary species" << std::
endl;
//
// Write other master species definitions, i.e. primary entities
//
chemistry_dat << "@species(H+, 1)" << std::endl;
for (cxxNameDouble::iterator iter = this->totals.begin();
iter != this->totals.end(); ++iter)
{
std::string name(iter->first);
if (name == "H(1)" || name == "E" || name == "H" || name == "O"
|| name == "Charge")
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);
chemistry_dat << "@species(" << s_ptr->name << ", " << s_ptr->
z << ")" << std::endl;
// regular mole balance
chemistry_dat << "@primary_entity(" << s_ptr->
name << ", 1e-9, liter, 1.0e-9)" << std::endl;
}
chemistry_dat << std::endl;
}
void
cxxSystem::ORCH_write_chemistry_total_O_H(std::ostream & chemistry_dat)
{
chemistry_dat << std::
endl << "//********* Define total hydrogen and oxygen" << std::endl;
// Write H total equation
chemistry_dat << "@var: total_hydrogen 0" << std::endl;
chemistry_dat << "@calc: (5, \"total_hydrogen = 0";
int i;
for (i = 0; i < count_s_x; i++)
{
// write in terms of orchestra components
if (s_x[i]->primary != NULL
|| (s_x[i]->secondary != NULL
&& s_x[i]->secondary->in == TRUE))
{
if (s_x[i]->h != 0)
{
chemistry_dat << "+";
if (s_x[i]->h != 1)
{
chemistry_dat << s_x[i]->h << "*";
}
chemistry_dat << "{" << s_x[i]->name << ".liter}";
}
}
}
chemistry_dat << "\")" << std::endl;
// Write O total equation
chemistry_dat << "@var: total_oxygen 0" << std::endl;
chemistry_dat << "@calc: (5, \"total_oxygen = 0";
for (i = 0; i < count_s_x; i++)
{
if (s_x[i]->o != 0)
{
// write in terms of orchestra components
if (s_x[i]->primary != NULL
|| (s_x[i]->secondary != NULL
&& s_x[i]->secondary->in == TRUE))
{
chemistry_dat << "+";
if (s_x[i]->o != 1)
{
chemistry_dat << s_x[i]->o << "*";
}
chemistry_dat << "{" << s_x[i]->name << ".liter}";
}
}
}
chemistry_dat << "\")" << std::endl;
chemistry_dat << std::endl;
}
void
cxxSystem::ORCH_write_input(std::ostream & input_dat)
{
//
// Write orchestra input file info
//
std::ostringstream headings, data;
data.precision(DBL_DIG - 1);
headings << "var: ";
data << "data: ";
headings << "tempc\t";
data << this->solution->get_tc() << "\t";
headings << "pH\t";
data << this->solution->get_ph() << "\t";
headings << "pe\t";
data << this->solution->get_pe() << "\t";
headings << "H2O.act\t";
data << 1 << "\t";
headings << "I\t";
data << this->solution->get_mu() << "\t";
for (cxxNameDouble::iterator iter = this->orch_totals.begin();
iter != this->orch_totals.end(); ++iter)
{
headings << iter->first << ".liter" << "\t";
data << iter->second << "\t";
}
// activity estimate
for (cxxNameDouble::iterator iter = this->totals.begin();
iter != this->totals.end(); ++iter)
{
if (iter->first == "O" || iter->first == "Charge")
continue;
std::string master_name;
struct master *master_ptr;
master_ptr = master_bsearch(iter->first);
assert(master_ptr != NULL);
cxxNameDouble ma = this->solution->get_master_activity();
cxxNameDouble::iterator it = ma.find(iter->first);
if (it == ma.end())
{
it = ma.find(master_ptr->s->secondary->elt->name);
}
headings << master_ptr->s->name << ".act\t";
if (it != ma.end())
{
data << pow(10., it->second) << "\t";
}
else
{
data << 1e-9 << "\t";
}
}
// Isotopes
//s_oss << "-Isotopes" << std::endl;
/*
{
for (std::list<cxxSolutionIsotope>::const_iterator it = this->isotopes.begin(); it != isotopes.end(); ++it) {
it->dump_raw(s_oss, indent + 2);
}
}
*/
// Write data to string
input_dat << headings.str() << std::endl;
input_dat << data.str() << std::endl;
return;
}
void
cxxSystem::ORCH_components(void)
{
// translate from H, O, charge to H+tot, e-tot, and H2Otot
cxxNameDouble::iterator it;
cxxNameDouble temp_totals;
//
// Set names in orch_totals to master species names
//
this->orch_totals.clear();
temp_totals = this->totals;
temp_totals.erase("H");
temp_totals.erase("O");
temp_totals.erase("Charge");
//temp_totals.dump_raw(std::cerr, 1);
for (it = temp_totals.begin(); it != temp_totals.end(); it++)
{
struct element *elt_ptr;
elt_ptr = element_store(it->first);
assert(elt_ptr != NULL);
struct master *master_ptr;
master_ptr = elt_ptr->primary;
assert(master_ptr != NULL);
double coef = master_ptr->coef;
if (master_ptr->coef == 0)
{
coef = 1;
}
this->orch_totals[master_ptr->s->name] = it->second / coef;
}
//
// Calculate h2otot
//
it = this->totals.find("O");
assert(it != this->totals.end());
double h2otot = it->second;
// subtract O in master species
for (it = temp_totals.begin(); it != temp_totals.end(); it++)
{
struct element *elt_ptr;
elt_ptr = element_store(it->first);
struct master *master_ptr;
master_ptr = elt_ptr->primary;
h2otot -= it->second * master_ptr->s->o;
}
//
// Calculate htot
//
it = this->totals.find("H");
assert(it != this->totals.end());
double htot = it->second - 2 * h2otot;
// subtract O in master species
for (it = temp_totals.begin(); it != temp_totals.end(); it++)
{
struct element *elt_ptr;
elt_ptr = element_store(it->first);
struct master *master_ptr;
master_ptr = elt_ptr->primary;
htot -= it->second * master_ptr->s->h;
}
//
// Calculate etot
//
it = this->totals.find("Charge");
assert(it != this->totals.end());
double etot = -it->second + htot;
// subtract O in master species
for (it = temp_totals.begin(); it != temp_totals.end(); it++)
{
struct element *elt_ptr;
elt_ptr = element_store(it->first);
struct master *master_ptr;
master_ptr = elt_ptr->primary;
etot += it->second * master_ptr->s->z;
}
//
// store h2otot, htot, etot in orch_totals
//
this->orch_totals["H2O"] = h2otot;
this->orch_totals["H+"] = htot;
this->orch_totals["e-"] = etot;
this->orch_totals.dump_raw(std::cerr, 1);
}
void
cxxSystem::ORCH_write_output_vars(std::ostream & outstream)
{
outstream << "\tnr_iter";
//
// 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.con";
//total_alkalinity;
outstream << "\tAlkalinity";
//orchestra master variables
outstream << "\tH+.diss";
outstream << "\te-.diss";
outstream << "\tH2O.diss";
//
// Write totals
for (cxxNameDouble::iterator it = this->orch_totals.begin();
it != this->orch_totals.end(); it++)
{
if (it->first == "H+" || it->first == "H2O" || it->first == "e-")
continue;
outstream << "\t" << it->first << ".diss";
}
outstream << "\tend_totals";
for (cxxNameDouble::iterator it = this->orch_totals.begin();
it != this->orch_totals.end(); it++)
{
if (it->first == "H+" || it->first == "H2O" || it->first == "e-")
continue;
outstream << "\t" << it->first << ".act";
}
outstream << "\tend_master_activities";
//
// Write all species activities and concentrations
//
int i;
for (i = 0; i < count_s_x; i++)
{
std::string name(s_x[i]->name);
std::replace(name.begin(), name.end(), '(', '[');
std::replace(name.begin(), name.end(), ')', ']');
outstream << "\t" << name.c_str() << ".act" << "\t" << name.
c_str() << ".con";
}
outstream << "\tend_species";
}
#endif /* */

13
System.h
View File

@ -12,6 +12,7 @@
//#include "Temperature.h"
#include "NameDouble.h"
#include "PHRQ_base.h"
#include "phreeqc_class.h"
class cxxSolution;
class cxxExchange;
class cxxGasPhase;
@ -75,17 +76,6 @@ public:
return this->totals;
}
#ifdef ORCHESTRA
void ORCH_components();
void ORCH_write(std::ostream & chemistry_dat, std::ostream & input_dat,
std::ostream & output_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_output_vars(std::ostream & outstream);
void ORCH_write_input(std::ostream & input_dat);
#endif /* */
private:
cxxSolution * solution;
cxxExchange * exchange;
@ -98,7 +88,6 @@ private:
cxxReaction * reaction;
cxxTemperature * temperature;
cxxNameDouble totals;
cxxNameDouble orch_totals;
};