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

Works with reaction calculation of pure phase equilibrium.

Browse Source 
Worked out relation between phreeqc components H,O,charge and Orchestra H+, e-, H2O.

It's a lot of code and still have not implemented SURFACE, EXCHANGE, GAS_PHASE,

much les KINETICS.

Convergence problem when including pyrite.

Also have REACTION (REQCTION_TEMP?);

git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@2181 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2007-08-16 18:51:35 +00:00
parent c78dd87daf
commit 7ab1e21ce7
13 changed files with 928 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Exchange.cxx
View File

@ -406,6 +406,7 @@ void cxxExchange::totalize()
for (std::list<cxxExchComp>::const_iterator it = exchComps.begin(); it != exchComps.end(); ++it)
{
this->totals.add_extensive(it->get_totals(), 1.0);
this->totals.add("Charge", it->get_charge_balance());
}
return;
}

70
ISolution.cxx
View File

@ -16,7 +16,7 @@
#include <cassert> // assert
#include <algorithm> // std::sort
#include <sstream>
extern void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
@ -426,7 +426,7 @@ void cxxISolution::ORCH_write_chemistry(std::ostream &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);
this->ORCH_write_chemistry_species(chemistry_dat);
ORCH_write_chemistry_species(chemistry_dat);
this->ORCH_write_chemistry_minerals(chemistry_dat);
}
@ -587,7 +587,6 @@ void cxxISolution::ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat)
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;
@ -656,41 +655,7 @@ void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat)
}
}
}
void cxxISolution::ORCH_write_chemistry_species(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Aqueous species" << std::endl;
//
// Write species definitions
//
for (int i = 0; i < count_s_x; i++) {
if (s_x[i]->primary != NULL) continue;
if (s_x[i]->secondary != NULL && s_x[i]->secondary->in == TRUE) continue;
// write species
std::string name(s_x[i]->name);
std::replace(name.begin(), name.end(), '(', '[');
std::replace(name.begin(), name.end(), ')', ']');
chemistry_dat << "@species(" << name.c_str() << ", " << s_x[i]->z << ")" << std::endl;
// write reaction
chemistry_dat << "@reaction(" << name.c_str() << ", " << pow(10.0, s_x[i]->lk);
struct rxn_token *next_token = s_x[i]->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 << std::endl;
}
void cxxISolution::ORCH_write_chemistry_minerals(std::ostream &chemistry_dat)
{
chemistry_dat << std::endl << "//********* Adjustments to mineral equilibrium" << std::endl;
@ -881,12 +846,13 @@ void cxxISolution::ORCH_write_output_vars(std::ostream &outstream)
//cb
outstream << "\tchargebalance";
//mass_water;
outstream << "\tH2O.diss";
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)
{
@ -930,3 +896,31 @@ void cxxISolution::ORCH_write_output_vars(std::ostream &outstream)
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();
}

2
ISolution.h
View File

@ -41,6 +41,7 @@ public:
//void dump_xml(std::ostream& os, unsigned int indent = 0)const;
void ConvertUnits();
void ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat);
void ORCH_write_chemistry(std::ostream &chemistry_dat);
void ORCH_write_input(std::ostream &input_dat);
void ORCH_write_output_vars(std::ostream &input_dat);
@ -61,7 +62,6 @@ private:
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);
};

72
PPassemblage.cxx
View File

@ -285,4 +285,76 @@ void cxxPPassemblage::add(const cxxPPassemblage &addee, double extensive)
//cxxNameDouble eltList;
this->eltList.add_extensive(addee.eltList, extensive);
}
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;
}
}
}

14
PPassemblage.h
View File

@ -31,23 +31,27 @@ public:
void read_raw(CParser& parser);
const cxxNameDouble& get_totals()const
{
return this->totals;
};
#ifdef USE_MPI
void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles);
void mpi_unpack(int *ints, int *ii, double *doubles, int *dd);
#endif
void totalize();
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);
private:
void add(const cxxPPassemblage &addee, double extensive);
// not written
void dump_xml(std::ostream& os, unsigned int indent = 0)const;
const cxxNameDouble& get_totals()const
{
return this->totals;
};
protected:
std::list<cxxPPassemblageComp> ppAssemblageComps;
cxxNameDouble eltList;

2
Solution.cxx
View File

@ -976,9 +976,11 @@ void cxxSolution::ORCH_read(std::vector <std::pair <std::string, double> > outpu
this->total_o = it->second; it++;
this->cb = it->second; it++;
this->mass_water = it->second * gfw_water; it++;
this->mass_water = 1.0;
this->total_alkalinity = it->second; it++;
it++; //orch total H+
it++; //orch total e-
it++; //orch total H2O
//cxxNameDouble totals;
char token[MAX_LENGTH];
while (it->first.compare("end_totals") != 0)

3
Solution.h
View File

@ -62,10 +62,11 @@ public:
double get_total_element(char *string)const;
void set_total(char *string, double value);
const cxxNameDouble& get_totals(void)const {return this->totals;}
void set_totals(cxxNameDouble &nd) {this->totals = nd; this->totals.type = cxxNameDouble::ND_ELT_MOLES;}
void clear_totals() {this->totals.clear();}
const cxxNameDouble& get_master_activity(void)const {return this->master_activity;}
double get_master_activity(char *string)const;
void set_master_activity(char *string, double value);

172
StorageBin.cxx
View File

@ -9,7 +9,7 @@
#include <iostream> // std::cout std::cerr
#include "Utils.h" // define first
#include "StorageBin.h"
#include "Solution.h"
#include "System.h"
#define EXTERNAL extern
#include "global.h"
#include "phqalloc.h"
@ -25,7 +25,73 @@ cxxStorageBin::cxxStorageBin()
{
// default constructor for cxxStorageBin
}
cxxStorageBin::cxxStorageBin(struct Use *use_ptr)
{
//Construct from use pointer
// Solution
if (use_ptr->solution_ptr != NULL)
{
cxxSolution entity(use_ptr->solution_ptr);
this->setSolution(use_ptr->n_solution_user, &entity);
}
// Exchange
if (use_ptr->exchange_ptr != NULL)
{
cxxExchange entity(use_ptr->exchange_ptr);
this->setExchange(use_ptr->n_exchange_user, &entity);
}
// gas_phase
if (use_ptr->gas_phase_ptr != NULL)
{
cxxGasPhase entity(use_ptr->gas_phase_ptr);
this->setGasPhase(use_ptr->n_gas_phase_user, &entity);
}
// kinetics
if (use_ptr->kinetics_ptr != NULL)
{
cxxKinetics entity(use_ptr->kinetics_ptr);
this->setKinetics(use_ptr->n_kinetics_user, &entity);
}
// pp_assemblage
if (use_ptr->pp_assemblage_ptr != NULL)
{
cxxPPassemblage entity(use_ptr->pp_assemblage_ptr);
this->setPPassemblage(use_ptr->n_pp_assemblage_user, &entity);
}
// s_s_assemblage
if (use_ptr->s_s_assemblage_ptr != NULL)
{
cxxSSassemblage entity(use_ptr->s_s_assemblage_ptr);
this->setSSassemblage(use_ptr->n_s_s_assemblage_user, &entity);
}
// surface
if (use_ptr->surface_ptr != NULL)
{
cxxSurface entity(use_ptr->surface_ptr);
this->setSurface(use_ptr->n_surface_user, &entity);
}
// mix
if (use_ptr->mix_ptr != NULL)
{
cxxMix entity(use_ptr->mix_ptr);
this->setMix(use_ptr->n_mix_user, &entity);
}
// reaction
if (use_ptr->irrev_ptr != NULL)
{
cxxReaction entity(use_ptr->irrev_ptr);
this->setReaction(use_ptr->n_irrev_user, &entity);
}
// reaction temperature
if (use_ptr->temperature_ptr != NULL)
{
cxxTemperature entity(use_ptr->temperature_ptr);
this->setTemperature(use_ptr->n_temperature_user, &entity);
}
// set system
this->setSystem(use_ptr);
}
cxxStorageBin::~cxxStorageBin()
{
}
@ -917,3 +983,107 @@ cxxExchange *cxxStorageBin::mix_cxxExchange(cxxMix &mixmap)
return(new_exch_ptr);
}
#endif
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);
}
void cxxStorageBin::setSystem(struct Use *use_ptr)
{
// Initialize
this->system.Initialize();
// Solution
if (use_ptr->solution_ptr != NULL)
{
std::map <int, cxxSolution>::iterator it = this->Solutions.find(use_ptr->n_solution_user);
if (it != this->Solutions.end())
{
this->system.setSolution(&(it->second));
}
}
// Exchange
if (use_ptr->exchange_ptr != NULL)
{
std::map <int, cxxExchange>::iterator it = this->Exchangers.find(use_ptr->n_exchange_user);
if (it != this->Exchangers.end())
{
this->system.setExchange(&(it->second));
}
}
// gas_phase
if (use_ptr->gas_phase_ptr != NULL)
{
std::map <int, cxxGasPhase>::iterator it = this->GasPhases.find(use_ptr->n_gas_phase_user);
if (it != this->GasPhases.end())
{
this->system.setGasPhase(&(it->second));
}
}
// kinetics
if (use_ptr->kinetics_ptr != NULL)
{
std::map <int, cxxKinetics>::iterator it = this->Kinetics.find(use_ptr->n_kinetics_user);
if (it != this->Kinetics.end())
{
this->system.setKinetics(&(it->second));
}
}
// pp_assemblage
if (use_ptr->pp_assemblage_ptr != NULL)
{
std::map <int, cxxPPassemblage>::iterator it = this->PPassemblages.find(use_ptr->n_pp_assemblage_user);
if (it != this->PPassemblages.end())
{
this->system.setPPassemblage(&(it->second));
}
}
// s_s_assemblage
if (use_ptr->s_s_assemblage_ptr != NULL)
{
std::map <int, cxxSSassemblage>::iterator it = this->SSassemblages.find(use_ptr->n_s_s_assemblage_user);
if (it != this->SSassemblages.end())
{
this->system.setSSassemblage(&(it->second));
}
}
// surface
if (use_ptr->surface_ptr != NULL)
{
std::map <int, cxxSurface>::iterator it = this->Surfaces.find(use_ptr->n_surface_user);
if (it != this->Surfaces.end())
{
this->system.setSurface(&(it->second));
}
}
// mix
if (use_ptr->mix_ptr != NULL)
{
std::map <int, cxxMix>::iterator it = this->Mixes.find(use_ptr->n_mix_user);
if (it != this->Mixes.end())
{
this->system.setMix(&(it->second));
}
}
// reaction
if (use_ptr->irrev_ptr != NULL)
{
std::map <int, cxxReaction>::iterator it = this->Reactions.find(use_ptr->n_irrev_user);
if (it != this->Reactions.end())
{
this->system.setReaction(&(it->second));
}
}
// reaction temperature
if (use_ptr->temperature_ptr != NULL)
{
std::map <int, cxxTemperature>::iterator it = this->Temperatures.find(use_ptr->n_temperature_user);
if (it != this->Temperatures.end())
{
this->system.setTemperature(&(it->second));
}
}
}

53
StorageBin.h
View File

@ -10,6 +10,7 @@
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "Surface.h"
#include "System.h"
#include "cxxMix.h"
#include "Reaction.h"
#include "Temperature.h"
@ -25,7 +26,7 @@ class cxxStorageBin
public:
cxxStorageBin();
cxxStorageBin(struct Use *use_ptr);
~cxxStorageBin();
void import_phreeqc(void);
@ -136,7 +137,51 @@ public:
Surfaces.erase(n_user);
}
void dump_raw(std::ostream& s_oss, unsigned int indent)const;
cxxMix *getMix(int n_user) {
if (this->Mixes.find(n_user) != this->Mixes.end()) {
return(&(this->Mixes.find(n_user)->second));
}
return (NULL);
}
void setMix(int n_user, cxxMix *entity) {
if (entity == NULL) return;
Mixes[n_user] = *entity;
}
void removeMix(int n_user) {
Mixes.erase(n_user);
}
cxxReaction *getReaction(int n_user) {
if (this->Reactions.find(n_user) != this->Reactions.end()) {
return(&(this->Reactions.find(n_user)->second));
}
return (NULL);
}
void setReaction(int n_user, cxxReaction *entity) {
if (entity == NULL) return;
Reactions[n_user] = *entity;
}
void removeReaction(int n_user) {
Reactions.erase(n_user);
}
cxxTemperature *getTemperature(int n_user) {
if (this->Temperatures.find(n_user) != this->Temperatures.end()) {
return(&(this->Temperatures.find(n_user)->second));
}
return (NULL);
}
void setTemperature(int n_user, cxxTemperature *entity) {
if (entity == NULL) return;
Temperatures[n_user] = *entity;
}
void removeTemperature(int n_user) {
Temperatures.erase(n_user);
}
void setSystem(struct Use *use_ptr);
void dump_raw(std::ostream& s_oss, unsigned int indent)const;
void dump_raw(std::ostream& s_oss, int i, unsigned int indent);
@ -161,7 +206,7 @@ public:
void mpi_send(int n, int task_number);
void mpi_recv(int task_number);
#endif
void ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat);
protected:
// Tidied classes
std::map<int, cxxSolution> Solutions;
@ -180,7 +225,7 @@ protected:
std::map<int, cxxMix> Mixes;
std::map<int, cxxReaction> Reactions;
std::map<int, cxxTemperature> Temperatures;
cxxSystem system;
public:
//static std::map<int, cxxStorageBin>& map;

1
Surface.cxx
View File

@ -680,6 +680,7 @@ void cxxSurface::totalize()
for (std::list<cxxSurfaceComp>::const_iterator it = surfaceComps.begin(); it != surfaceComps.end(); ++it)
{
this->totals.add_extensive(it->get_totals(), 1.0);
this->totals.add("Charge", it->get_charge_balance());
}
return;
}

2
SurfaceComp.h
View File

@ -27,6 +27,8 @@ public:
char *get_phase_name()const {return this->phase_name;}
char *get_rate_name()const {return this->rate_name;}
char *get_formula()const {return this->formula;}
double get_charge_balance()const { return this->charge_balance;}
void set_charge_balance(double d) { this->charge_balance = d;}
static struct surface_comp *cxxSurfaceComp2surface_comp(std::list<cxxSurfaceComp>& el);

513
System.cxx Normal file
View File

@ -0,0 +1,513 @@
#include "System.h"
#include <algorithm> // std::replace
extern void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
cxxSystem::cxxSystem(void)
{
this->solution = NULL;
this->exchange = NULL;
this->ppassemblage = NULL;
this->gasphase = NULL;
this->ssassemblage = NULL;
this->kinetics = NULL;
this->surface = NULL;
this->mix = NULL;
this->reaction = NULL;
this->temperature = NULL;
}
cxxSystem::~cxxSystem(void)
{
}
void cxxSystem::Initialize(void)
{
this->solution = NULL;
this->exchange = NULL;
this->ppassemblage = NULL;
this->gasphase = NULL;
this->ssassemblage = NULL;
this->kinetics = NULL;
this->surface = NULL;
this->mix = NULL;
this->reaction = NULL;
this->temperature = NULL;
}
void cxxSystem::totalize(void)
{
//initialize
this->totals.clear();
//add solution
if (this->solution != NULL)
{
this->totals["O"] = this->solution->get_total_o();
this->totals["H"] = this->solution->get_total_h();
this->totals["Charge"] = this->solution->get_cb();
this->totals.add_extensive(this->solution->get_totals(), 1.0);
}
if (this->exchange != NULL)
{
this->exchange->totalize();
this->totals.add_extensive(this->exchange->get_totals(), 1.0);
}
if (this->ppassemblage != NULL)
{
this->ppassemblage->totalize();
this->totals.add_extensive(this->ppassemblage->get_totals(), 1.0);
}
if (this->gasphase != NULL)
{
this->gasphase->totalize();
this->totals.add_extensive(this->gasphase->get_totals(), 1.0);
}
if (this->ssassemblage != NULL)
{
this->ssassemblage->totalize();
this->totals.add_extensive(this->ssassemblage->get_totals(), 1.0);
}
if (this->surface != NULL)
{
this->ssassemblage->totalize();
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;
}
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;
}
#ifdef SKIP
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";
for (cxxNameDouble::iterator iter = this->totals.begin(); iter != this->totals.end(); ++iter)
{
if (iter->first == "O") continue;
std::string master_name;
struct master *master_ptr;
master_ptr = master_bsearch (iter->first);
assert (master_ptr != NULL);
double coef = master_ptr->coef;
if (master_ptr->coef == 0)
{
coef = 1;
}
if (iter->first == "H")
{
headings << "system_hydrogen" << "\t";
data << iter->second << "\t";
continue;
} else
{
headings << master_ptr->s->name << ".liter" << "\t";
data << iter->second / coef << "\t";
}
// activity estimate
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;
}
#endif
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;
struct species *s_ptr = s_search(it->first);
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";
}

73
System.h Normal file
View File

@ -0,0 +1,73 @@
#if !defined(SYSTEM_H_INCLUDED)
#define SYSTEM_H_INCLUDED
#include "Solution.h"
#include "Exchange.h"
#include "GasPhase.h"
#include "cxxKinetics.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "Surface.h"
#include "cxxMix.h"
#include "Reaction.h"
#include "Temperature.h"
class cxxSystem
{
public:
cxxSystem(void);
public:
~cxxSystem(void);
void Initialize(void);
void setSolution(cxxSolution *entity) {
this->solution = entity;
}
void setExchange(cxxExchange *entity) {
this->exchange = entity;
}
void setPPassemblage(cxxPPassemblage *entity) {
this->ppassemblage = entity;
}
void setGasPhase(cxxGasPhase *entity) {
this->gasphase = entity;
}
void setSSassemblage(cxxSSassemblage *entity) {
this->ssassemblage = entity;
}
void setKinetics(cxxKinetics *entity) {
this->kinetics = entity;
}
void setSurface(cxxSurface *entity) {
this->surface = entity;
}
void setMix(cxxMix *entity) {
this->mix = entity;
}
void setReaction(cxxReaction *entity) {
this->reaction = entity;
}
void setTemperature(cxxTemperature *entity) {
this->temperature = entity;
}
void totalize();
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);
private:
cxxSolution *solution;
cxxExchange *exchange;
cxxPPassemblage *ppassemblage;
cxxGasPhase *gasphase;
cxxSSassemblage *ssassemblage;
cxxKinetics *kinetics;
cxxSurface *surface;
cxxMix *mix;
cxxReaction *reaction;
cxxTemperature *temperature;
cxxNameDouble totals;
cxxNameDouble orch_totals;
};
#endif // !defined(SYSTEM_H_INCLUDED)