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

Beginning to revise cpp classes.

Browse Source 
Have worked on mixing solutions.

Roughed in mixing Exchange.h



git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@1722 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2007-01-30 01:05:44 +00:00
parent bd852de57c
commit 93a7759e1b
9 changed files with 286 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
ExchComp.cxx
View File

@ -102,7 +102,7 @@ cxxExchComp::cxxExchComp(std::vector<cxxExchComp> &ec_vector, std::vector<double
this->la += it_ec->la*intensive; this->la += it_ec->la*intensive;
this->charge_balance += it_ec->charge_balance*extensive; this->charge_balance += it_ec->charge_balance*extensive;
this->phase_proportion += it_ec->phase_proportion*intensive; this->phase_proportion += it_ec->phase_proportion*intensive;
this->totals.add(it_ec->totals, extensive); this->totals.add_extensive(it_ec->totals, extensive);
it_ec++; it_ec++;
it_f++; it_f++;
} }

14
Exchange.h
View File

@ -1,17 +1,15 @@
#if !defined(EXCHANGE_H_INCLUDED) #if !defined(EXCHANGE_H_INCLUDED)
#define EXCHANGE_H_INCLUDED #define EXCHANGE_H_INCLUDED
#include "NumKeyword.h"
#define EXTERNAL extern #define EXTERNAL extern
#include "global.h" #include "global.h"
#include <cassert> // assert #include <cassert> // assert
#include <map> // std::map #include <map> // std::map
#include <string> // std::string #include <string> // std::string
#include <list> // std::list #include <list> // std::list
#include <vector> // std::vector
#include "char_star.h"
#include "ExchComp.h" #include "ExchComp.h"
#include "cxxMix.h"
class cxxExchange : public cxxNumKeyword class cxxExchange : public cxxNumKeyword
{ {
@ -19,14 +17,13 @@ class cxxExchange : public cxxNumKeyword
public: public:
cxxExchange(); cxxExchange();
cxxExchange(struct exchange *); cxxExchange(struct exchange *);
cxxExchange(const std::map<int, cxxExchange> &exchange_map, cxxMix &mx);
~cxxExchange(); ~cxxExchange();
struct exchange *cxxExchange2exchange(); struct exchange *cxxExchange2exchange();
struct exch_comp *cxxExchComp2exch_comp(); struct exch_comp *cxxExchComp2exch_comp();
void dump_xml(std::ostream& os, unsigned int indent = 0)const;
void dump_raw(std::ostream& s_oss, unsigned int indent)const; void dump_raw(std::ostream& s_oss, unsigned int indent)const;
void read_raw(CParser& parser); void read_raw(CParser& parser);
@ -50,6 +47,11 @@ public:
void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles); void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles);
void mpi_unpack(int *ints, int *ii, double *doubles, int *dd); void mpi_unpack(int *ints, int *ii, double *doubles, int *dd);
#endif #endif
private:
void zero();
void add(const cxxExchange &ex, double intensive, double extensive);
// not written
void dump_xml(std::ostream& os, unsigned int indent = 0)const;
protected: protected:
std::list<cxxExchComp> exchComps; std::list<cxxExchComp> exchComps;

68
NameDouble.cxx
View File

@ -275,21 +275,63 @@ CParser::STATUS_TYPE cxxNameDouble::read_raw(CParser& parser, std::istream::pos_
return CParser::PARSER_OK; return CParser::PARSER_OK;
} }
void cxxNameDouble::add(const cxxNameDouble &old, double factor) void cxxNameDouble::add_extensive(const cxxNameDouble &addee, double factor)
// //
// constructor for cxxNameDouble from list of elt_list // Sums two name doubles, this + factor*nd2
// //
{ {
for (cxxNameDouble::const_iterator it = old.begin(); it != old.end(); it++) { if (factor == 0) return;
cxxNameDouble::iterator current = (*this).find(it->first); assert (factor > 0.0);
if (current != (*this).end()) { for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++)
(*this)[it->first] = current->second + it->second * factor; {
} else { cxxNameDouble::iterator current = (*this).find(it->first);
(*this)[it->first] = it->second * factor; if (current != (*this).end())
} {
} (*this)[it->first] = current->second + it->second * factor;
} else {
(*this)[it->first] = it->second * factor;
}
}
} }
void cxxNameDouble::add_intensive(const cxxNameDouble &addee, double f1, double f2)
//
// Sums two name doubles, this*f1 + f2*nd2
//
{
assert(f1 > 0 && f2 > 0);
for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++)
{
cxxNameDouble::iterator current = (*this).find(it->first);
if (current != (*this).end())
{
(*this)[it->first] = f1*current->second + f2*it->second;
} else {
(*this)[it->first] = f2 * it->second;
}
}
}
void cxxNameDouble::add_log_activities(const cxxNameDouble &addee, double f1, double f2)
//
// Sums two name doubles, this*f1 + f2*nd2, assuming log values
//
{
assert (f1 > 0 && f2 > 0);
for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end(); it++)
{
cxxNameDouble::iterator current = (*this).find(it->first);
if (current != (*this).end())
{
double a1 = pow(10., current->second);
double a2 = pow(10., it->second);
(*this)[it->first] = log10(f1*a1 + f2*a2);
} else {
//double a2 = pow(10. it->second);
//(*this)[it->first] = log10(f2 * a2);
(*this)[it->first] = it->second + log10(f2);
}
}
}
void cxxNameDouble::add(char * key, double total) void cxxNameDouble::add(char * key, double total)
// //
// add to total for a specified element // add to total for a specified element

4
NameDouble.h
View File

@ -44,7 +44,9 @@ public:
CParser::STATUS_TYPE read_raw(CParser& parser, std::istream::pos_type& pos); CParser::STATUS_TYPE read_raw(CParser& parser, std::istream::pos_type& pos);
void add(const cxxNameDouble &old, double factor); void add_extensive(const cxxNameDouble &old, double factor);
void add_intensive(const cxxNameDouble &addee, double fthis, double f2);
void add_log_activities(const cxxNameDouble &addee, double fthis, double f2);
void add(char * key, double total); void add(char * key, double total);
void insert(char *str, double d) { void insert(char *str, double d) {

172
Solution.cxx
View File

@ -27,20 +27,21 @@ cxxSolution::cxxSolution()
// //
: cxxNumKeyword() : cxxNumKeyword()
{ {
tc = 25.0; this->tc = 25.0;
ph = 7.0; this->ph = 7.0;
pe = 4.0; this->pe = 4.0;
mu = 1e-7; this->mu = 1e-7;
ah2o = 1.0; this->ah2o = 1.0;
total_h = 111.1; this->total_h = 111.1;
total_o = 55.55; this->total_o = 55.55;
cb = 0.0; this->cb = 0.0;
mass_water = 1.0; this->mass_water = 1.0;
total_alkalinity = 0.0; this->total_alkalinity = 0.0;
totals.type = cxxNameDouble::ND_ELT_MOLES; this->totals.type = cxxNameDouble::ND_ELT_MOLES;
master_activity.type = cxxNameDouble::ND_SPECIES_LA; this->master_activity.type = cxxNameDouble::ND_SPECIES_LA;
species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; this->species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA;
} }
/*
cxxSolution::cxxSolution(double zero) cxxSolution::cxxSolution(double zero)
// //
// empty cxxSolution constructor // empty cxxSolution constructor
@ -62,6 +63,7 @@ cxxSolution::cxxSolution(double zero)
master_activity.type = cxxNameDouble::ND_SPECIES_LA; master_activity.type = cxxNameDouble::ND_SPECIES_LA;
species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA; species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA;
} }
*/
cxxSolution::cxxSolution(struct solution *solution_ptr) cxxSolution::cxxSolution(struct solution *solution_ptr)
// //
// constructor for cxxSolution from struct solution // constructor for cxxSolution from struct solution
@ -75,18 +77,18 @@ isotopes(solution_ptr)
{ {
this->set_description(solution_ptr->description); this->set_description(solution_ptr->description);
n_user = solution_ptr->n_user; this->n_user = solution_ptr->n_user;
n_user_end = solution_ptr->n_user_end; this->n_user_end = solution_ptr->n_user_end;
tc = solution_ptr->tc; this->tc = solution_ptr->tc;
ph = solution_ptr->ph; this->ph = solution_ptr->ph;
pe = solution_ptr->solution_pe; this->pe = solution_ptr->solution_pe;
mu = solution_ptr->mu; this->mu = solution_ptr->mu;
ah2o = solution_ptr->ah2o; this->ah2o = solution_ptr->ah2o;
total_h = solution_ptr->total_h; this->total_h = solution_ptr->total_h;
total_o = solution_ptr->total_o; this->total_o = solution_ptr->total_o;
cb = solution_ptr->cb; this->cb = solution_ptr->cb;
mass_water = solution_ptr->mass_water; this->mass_water = solution_ptr->mass_water;
total_alkalinity = solution_ptr->total_alkalinity; this->total_alkalinity = solution_ptr->total_alkalinity;
// Totals filled in constructor, just save description and moles // Totals filled in constructor, just save description and moles
@ -101,6 +103,83 @@ isotopes(solution_ptr)
// Species_gamma in constructor // Species_gamma in constructor
} }
#ifdef SKIP
cxxSolution::cxxSolution( const std::map<int, cxxSolution>& solutions, cxxMix &mix)
//
// constructor for cxxSolution from struct solution
//cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap)
{
/*
* mixes solutions based on cxxMix structure, returns new solution
* return solution must be freed by calling method
*/
double intensive, extensive;
/*
* Zero out global solution data
*/
this->zero();
/*
* Determine sum of mixing fractions
*/
extensive = 0.0;
std::map<int, double> *mixcomps = mix.comps();
std::map<int, double>::const_iterator it;
for (it = mixcomps->begin(); it != mixcomps->end(); it++) {
extensive += it->second;
}
/*
* Add solutions
*/
for (it = mixcomps->begin(); it != mixcomps->end(); it++)
{
const cxxSolution *cxxsoln_ptr1 = &(solutions.find(it->first)->second);
if (cxxsoln_ptr1 == NULL)
{
sprintf(error_string, "Solution %d not found in mix_cxxSolutions.", it->first);
error_msg(error_string, CONTINUE);
input_error++;
return;
}
intensive = it->second/extensive;
double ext = it->second;
this->add(*cxxsoln_ptr1, intensive, it->second);
}
}
#endif
cxxSolution::cxxSolution( const std::map<int, cxxSolution>& solutions, cxxMix &mix)
//
// constructor for cxxSolution from mixture of solutions
//
{
//
// Zero out solution data
//
this->zero();
//
// Mix solutions
//
std::map<int, double> *mixcomps = mix.comps();
std::map<int, double>::const_iterator it;
for (it = mixcomps->begin(); it != mixcomps->end(); it++)
{
const cxxSolution *cxxsoln_ptr1 = &(solutions.find(it->first)->second);
if (cxxsoln_ptr1 == NULL)
{
sprintf(error_string, "Solution %d not found in mix_cxxSolutions.", it->first);
error_msg(error_string, CONTINUE);
input_error++;
return;
}
this->add(*cxxsoln_ptr1, it->second);
}
}
#ifdef SKIP #ifdef SKIP
cxxSolution::cxxSolution(const cxxSolution &old, double intensive, double extensive) cxxSolution::cxxSolution(const cxxSolution &old, double intensive, double extensive)
// //
@ -669,6 +748,23 @@ void cxxSolution::read_raw(CParser& parser)
return; return;
} }
void cxxSolution::zero()
{
this->tc = 0.0;
this->ph = 0.0;
this->pe = 0.0;
this->mu = 0.0;
this->ah2o = 0.0;
this->total_h = 0.0;
this->total_o = 0.0;
this->cb = 0.0;
this->mass_water = 0.0;
this->total_alkalinity = 0.0;
this->totals.type = cxxNameDouble::ND_ELT_MOLES;
this->master_activity.type = cxxNameDouble::ND_SPECIES_LA;
this->species_gamma.type = cxxNameDouble::ND_SPECIES_GAMMA;
}
#ifdef SKIP
void cxxSolution::add(const cxxSolution &old, double intensive, double extensive) void cxxSolution::add(const cxxSolution &old, double intensive, double extensive)
// //
// Add existing solution to "this" solution // Add existing solution to "this" solution
@ -722,7 +818,31 @@ void cxxSolution::add(const cxxSolution &old, double intensive, double extensive
cxxNameDouble species_gamma; cxxNameDouble species_gamma;
*/ */
} }
#endif
void cxxSolution::add(const cxxSolution &addee, double extensive)
//
// Add existing solution to "this" solution
//
{
double ext1 = this->mass_water;
double ext2 = addee.mass_water * extensive;
double f1 = ext1/(ext1 + ext2);
double f2 = ext2/(ext1 + ext2);
this->tc = f1*this->tc + f2*addee.tc;
this->ph = f1*this->ph + f2*addee.ph;
this->pe = f1*this->pe + f2*addee.pe;
this->mu = f1*this->mu + f2*addee.mu;
this->ah2o = f1*this->mu + f2*addee.ah2o;
this->total_h += addee.total_h * extensive;
this->total_o += addee.total_o * extensive;
this->cb += addee.cb * extensive;
this->mass_water += addee.mass_water * extensive;
this->total_alkalinity += addee.total_alkalinity * extensive;
this->totals.add_extensive(addee.totals, extensive);
this->master_activity.add_log_activities(addee.master_activity, f1, f2);
this->species_gamma.add_intensive(addee.species_gamma, f1, f2);
this->isotopes.add(addee.isotopes, f2, extensive);
}
double cxxSolution::get_total(char *string)const double cxxSolution::get_total(char *string)const
{ {
cxxNameDouble::const_iterator it = this->totals.find(string); cxxNameDouble::const_iterator it = this->totals.find(string);

22
Solution.h
View File

@ -11,7 +11,7 @@
#include <cassert> // assert #include <cassert> // assert
#include <map> // std::map #include <map> // std::map
#include <string> // std::string #include <string> // std::string
#include <list> // std::list
#include <vector> // std::vector #include <vector> // std::vector
#include "char_star.h" #include "char_star.h"
@ -21,12 +21,11 @@ class cxxSolution : public cxxNumKeyword
public: public:
cxxSolution(); cxxSolution();
cxxSolution(double zero); //cxxSolution(double zero);
cxxSolution(struct solution *); cxxSolution(struct solution *);
cxxSolution(const std::map<int, cxxSolution> &solution_map, cxxMix &mx); cxxSolution(const std::map<int, cxxSolution> &solution_map, cxxMix &mx);
cxxSolution(const cxxSolution &old, double intensive, double extensive); //cxxSolution(const cxxSolution &old, double intensive, double extensive);
//cxxSolution(const cxxSolution&);
//cxxSolution(const cxxSolution&);
~cxxSolution(); ~cxxSolution();
//static cxxSolution& read(CParser& parser); //static cxxSolution& read(CParser& parser);
@ -77,13 +76,10 @@ public:
struct solution *cxxSolution2solution(); struct solution *cxxSolution2solution();
void dump_xml(std::ostream& os, unsigned int indent = 0)const; void dump_raw(std::ostream& s_oss, unsigned int indent)const;
void dump_raw(std::ostream& s_oss, unsigned int indent)const;
void read_raw(CParser& parser); void read_raw(CParser& parser);
void add(const cxxSolution &sol, double intensive, double extensive);
#ifdef USE_MPI #ifdef USE_MPI
void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles); void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles);
@ -91,6 +87,14 @@ public:
void mpi_send(int task_number); void mpi_send(int task_number);
void mpi_recv(int task_number); void mpi_recv(int task_number);
#endif #endif
private:
void zero();
//void add(const cxxSolution &sol, double intensive, double extensive);
void add(const cxxSolution &addee, double extensive);
// not checked
void dump_xml(std::ostream& os, unsigned int indent = 0)const;
protected: protected:
double tc; double tc;
double ph; double ph;

93
SolutionIsotopeList.cxx
View File

@ -32,54 +32,55 @@ cxxSolutionIsotopeList::cxxSolutionIsotopeList(struct solution *solution_ptr)
} }
void cxxSolutionIsotopeList::add(cxxSolutionIsotopeList old, double intensive, double extensive) void cxxSolutionIsotopeList::add(cxxSolutionIsotopeList old, double intensive, double extensive)
{ {
for (cxxSolutionIsotopeList::const_iterator itold = old.begin(); itold != old.end(); ++itold) {
for (cxxSolutionIsotopeList::const_iterator itold = old.begin(); itold != old.end(); ++itold) { //for (std::list <cxxSolutionIsotope>::const_iterator itold = old.isotopes.begin(); itold != old.isotopes.end(); ++itold) {
bool found = false;
//for (std::list <cxxSolutionIsotope>::const_iterator itold = old.isotopes.begin(); itold != old.isotopes.end(); ++itold) { for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it)
bool found = false; {
for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) { //for (std::list <cxxSolutionIsotope>::iterator it = this->isotopes.begin(); it != this->isotopes.end(); ++it) {
//for (std::list <cxxSolutionIsotope>::iterator it = this->isotopes.begin(); it != this->isotopes.end(); ++it) { if ((it->isotope_number == itold->isotope_number) &&
if ((it->isotope_number == itold->isotope_number) && (it->elt_name == itold->elt_name) &&
(it->elt_name == itold->elt_name) && (it->isotope_name == itold->isotope_name))
(it->isotope_name == itold->isotope_name)) { {
it->total += itold->total * extensive; it->total += itold->total * extensive;
it->ratio += itold->ratio * intensive; it->ratio += itold->ratio * intensive;
it->ratio_uncertainty += itold->ratio_uncertainty * intensive; it->ratio_uncertainty += itold->ratio_uncertainty * intensive;
it->ratio_uncertainty_defined = (it->ratio_uncertainty_defined || itold->ratio_uncertainty_defined); it->ratio_uncertainty_defined = (it->ratio_uncertainty_defined || itold->ratio_uncertainty_defined);
found = true; found = true;
break; break;
} }
} }
if (!found) { if (!found)
cxxSolutionIsotope iso; {
iso.total = itold->total * extensive; cxxSolutionIsotope iso;
iso.ratio = itold->ratio * intensive; iso.total = itold->total * extensive;
iso.ratio_uncertainty = itold->ratio_uncertainty * intensive; iso.ratio = itold->ratio * intensive;
iso.ratio_uncertainty_defined = itold->ratio_uncertainty_defined; iso.ratio_uncertainty = itold->ratio_uncertainty * intensive;
this->push_back(iso); iso.ratio_uncertainty_defined = itold->ratio_uncertainty_defined;
} this->push_back(iso);
} }
}
} }
struct isotope * cxxSolutionIsotopeList::cxxSolutionIsotopeList2isotope() struct isotope * cxxSolutionIsotopeList::cxxSolutionIsotopeList2isotope()
{ {
struct isotope *iso; struct isotope *iso;
if (this->size() <= 0) { if (this->size() <= 0) {
return NULL; return NULL;
} else { } else {
iso = (struct isotope *) PHRQ_malloc((size_t) ((this->size()) * sizeof(struct isotope))); iso = (struct isotope *) PHRQ_malloc((size_t) ((this->size()) * sizeof(struct isotope)));
if (iso == NULL) malloc_error(); if (iso == NULL) malloc_error();
int i = 0; int i = 0;
for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) { for (cxxSolutionIsotopeList::iterator it = this->begin(); it != this->end(); ++it) {
iso[i].isotope_number = it->isotope_number; iso[i].isotope_number = it->isotope_number;
iso[i].elt_name = it->elt_name; iso[i].elt_name = it->elt_name;
iso[i].total = it->total; iso[i].total = it->total;
iso[i].ratio = it->ratio; iso[i].ratio = it->ratio;
iso[i].ratio_uncertainty = it->ratio_uncertainty; iso[i].ratio_uncertainty = it->ratio_uncertainty;
iso[i].master = it->master(); iso[i].master = it->master();
iso[i].primary = it->primary(); iso[i].primary = it->primary();
i++; i++;
} }
} }
return(iso); return(iso);
} }

4
StorageBin.cxx
View File

@ -495,7 +495,7 @@ void cxxStorageBin::remove(int n)
// Surface // Surface
this->Surfaces.erase(n); this->Surfaces.erase(n);
} }
#ifdef SKIP
cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap) cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap)
{ {
@ -536,7 +536,7 @@ cxxSolution *cxxStorageBin::mix_cxxSolutions(cxxMix &mixmap)
} }
return(cxxsoln_ptr); return(cxxsoln_ptr);
} }
#endif
struct system *cxxStorageBin::cxxStorageBin2system(int n) struct system *cxxStorageBin::cxxStorageBin2system(int n)
// //
// make a system from storagebin // make a system from storagebin

14
StorageBin.h
View File

@ -42,8 +42,8 @@ public:
} }
return (NULL); return (NULL);
} }
void setSolution(int n_user, cxxSolution *entity) { void setSolution(int n_user, cxxSolution &entity) {
Solutions[n_user] = *entity; Solutions[n_user] = entity;
} }
void removeSolution(int n_user) { void removeSolution(int n_user) {
Solutions.erase(n_user); Solutions.erase(n_user);
@ -137,9 +137,17 @@ public:
struct system *cxxStorageBin2system(int i); struct system *cxxStorageBin2system(int i);
cxxSolution *mix_cxxSolutions(cxxMix &mixmap); //cxxSolution *mix_cxxSolutions(cxxMix &mixmap);
cxxExchange *mix_cxxExchange(cxxMix &mixmap); cxxExchange *mix_cxxExchange(cxxMix &mixmap);
const std::map<int, cxxSolution>& getSolutions()const {return this->Solutions;};
const std::map<int, cxxExchange>& getExchangers()const {return this->Exchangers;};
const std::map<int, cxxGasPhase>& getGasPhases()const {return this->GasPhases;};
const std::map<int, cxxKinetics>& getKinetics()const {return this->Kinetics;};
const std::map<int, cxxPPassemblage>& getPPassemblages()const {return this->PPassemblages;};
const std::map<int, cxxSSassemblage>& getSSassemblages()const {return this->SSassemblages;};
const std::map<int, cxxSurface>& getSurfaces()const {return this->Surfaces;};
#ifdef USE_MPI #ifdef USE_MPI
void mpi_send(int n, int task_number); void mpi_send(int n, int task_number);
void mpi_recv(int task_number); void mpi_recv(int task_number);