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

Merge commit 'b9f44da20a78a23395d38205c001298f6a68bd0c'

Browse Source
This commit is contained in:
Darth Vader 2021-03-05 23:49:46 +00:00
commit 996c35291b
8 changed files with 311 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

71
src/phreeqcpp/GasComp.cxx
View File

@ -30,6 +30,9 @@ cxxGasComp::cxxGasComp(PHRQ_io *io)
p_read = 0.0;
moles = 0.0;
initial_moles = 0.0;
p = 0.0;
phi = 0.0;
f = 0.0;
}
cxxGasComp::~cxxGasComp()
@ -56,6 +59,9 @@ cxxGasComp::dump_raw(std::ostream & s_oss, unsigned int indent) const
s_oss << indent0 << "# GasComp workspace variables #\n";
s_oss << indent0 << "-initial_moles " << this->initial_moles << "\n";
s_oss << indent0 << "-p " << this->p << "\n";
s_oss << indent0 << "-phi " << this->phi << "\n";
s_oss << indent0 << "-f " << this->f << "\n";
}
bool
@ -125,6 +131,37 @@ cxxGasComp::read_raw(CParser & parser, bool check)
PHRQ_io::OT_CONTINUE);
}
break;
case 5: // p
if (!(parser.get_iss() >> this->p))
{
this->p = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for pressure.",
PHRQ_io::OT_CONTINUE);
}
break;
case 6: // phi
if (!(parser.get_iss() >> this->phi))
{
this->phi = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for phi.",
PHRQ_io::OT_CONTINUE);
}
break;
case 7: // f
if (!(parser.get_iss() >> this->f))
{
this->f = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for f.",
PHRQ_io::OT_CONTINUE);
}
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
@ -150,9 +187,9 @@ cxxGasComp::add(const cxxGasComp & addee, LDBLE extensive)
if (addee.phase_name.size() == 0)
return;
/*
ext1 = this->moles;
ext2 = addee.moles * extensive;
double f1, f2;
double ext1 = this->moles;
double ext2 = addee.moles * extensive;
if (ext1 + ext2 != 0)
{
f1 = ext1 / (ext1 + ext2);
@ -163,13 +200,15 @@ cxxGasComp::add(const cxxGasComp & addee, LDBLE extensive)
f1 = 0.5;
f2 = 0.5;
}
*/
assert(this->phase_name == addee.phase_name);
this->p_read += addee.p_read * extensive;
this->p_read = this->p_read*f1 + addee.p_read * f2;
this->moles += addee.moles * extensive;
this->initial_moles += addee.initial_moles * extensive;
this->p = this->p * f1 + addee.p * f2;
this->phi = this->phi * f1 + addee.phi * f2;
this->f = this->f * f1 + addee.f * f2;
}
void
@ -178,6 +217,9 @@ cxxGasComp::multiply(LDBLE extensive)
this->p_read *= extensive;
this->moles *= extensive;
this->initial_moles *= extensive;
this->p *= 1.0;
this->phi *= 1.0;
this->f *= 1.0;
}
void
cxxGasComp::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles)
@ -186,6 +228,9 @@ cxxGasComp::Serialize(Dictionary & dictionary, std::vector < int >&ints, std::ve
doubles.push_back(this->moles);
doubles.push_back(this->p_read);
doubles.push_back(this->initial_moles);
doubles.push_back(this->p);
doubles.push_back(this->phi);
doubles.push_back(this->f);
}
void
@ -196,13 +241,19 @@ cxxGasComp::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->moles = doubles[dd++];
this->p_read = doubles[dd++];
this->initial_moles = doubles[dd++];
this->p = doubles[dd++];
this->phi = doubles[dd++];
this->f = doubles[dd++];
}
const std::vector< std::string >::value_type temp_vopts[] = {
std::vector< std::string >::value_type("phase_name"), // 0
std::vector< std::string >::value_type("name"), // 1
std::vector< std::string >::value_type("p_read"), // 2
std::vector< std::string >::value_type("moles"), // 3
std::vector< std::string >::value_type("initial_moles") // 4
std::vector< std::string >::value_type("phase_name"), // 0
std::vector< std::string >::value_type("name"), // 1
std::vector< std::string >::value_type("p_read"), // 2
std::vector< std::string >::value_type("moles"), // 3
std::vector< std::string >::value_type("initial_moles"), // 4
std::vector< std::string >::value_type("p"), // 5
std::vector< std::string >::value_type("phi"), // 6
std::vector< std::string >::value_type("f") // 7
};
const std::vector< std::string > cxxGasComp::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

9
src/phreeqcpp/GasComp.h
View File

@ -28,6 +28,12 @@ class cxxGasComp: public PHRQ_base
void Set_p_read(LDBLE t) {this->p_read = t;}
LDBLE Get_moles() const {return this->moles;}
void Set_moles(LDBLE t) {this->moles = t;}
LDBLE Get_p() const { return this->p; }
void Set_p(LDBLE t) { this->p = t; }
LDBLE Get_phi() const { return this->phi; }
void Set_phi(LDBLE t) { this->phi = t; }
LDBLE Get_f() const { return this->f; }
void Set_f(LDBLE t) { this->f = t; }
LDBLE Get_initial_moles() const {return this->initial_moles;}
void Set_initial_moles(LDBLE t) {this->initial_moles = t;}
@ -44,6 +50,9 @@ class cxxGasComp: public PHRQ_base
LDBLE p_read;
// internal workspace
LDBLE initial_moles;
LDBLE p;
LDBLE phi;
LDBLE f;
const static std::vector < std::string > vopts;
};

86
src/phreeqcpp/GasPhase.cxx
View File

@ -649,6 +649,92 @@ LDBLE cxxGasPhase::Calc_total_moles(void)const
}
return tot;
}
void cxxGasPhase::Delete_component(const std::string comp_name)
{
for (size_t i = 0; i < this->gas_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0)
{
this->gas_comps.erase(this->gas_comps.begin() + i); // To delete the ith element
break;
}
}
}
void cxxGasPhase::Set_component_moles(const std::string comp_name, const double moles)
{
if (moles < 0.0)
{
this->Delete_component(comp_name);
}
else
{
cxxGasComp* ptr = this->Find_comp(comp_name.c_str());
if (ptr != NULL)
{
ptr->Set_moles(moles);
}
else
{
cxxGasComp temp_comp;
temp_comp.Set_phase_name(comp_name);
temp_comp.Set_moles(moles);
this->gas_comps.push_back(temp_comp);
}
}
}
double cxxGasPhase::Get_component_moles(const std::string comp_name)
{
double moles = -1.0;
for (size_t i = 0; i < this->gas_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0)
{
moles = this->gas_comps[i].Get_moles();
break;
}
}
return moles;
}
double cxxGasPhase::Get_component_p(const std::string comp_name)
{
double p = -1.0;
for (size_t i = 0; i < this->gas_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0)
{
p = this->gas_comps[i].Get_p();
break;
}
}
return p;
}
double cxxGasPhase::Get_component_phi(const std::string comp_name)
{
double phi = -1.0;
for (size_t i = 0; i < this->gas_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0)
{
phi = this->gas_comps[i].Get_phi();
break;
}
}
return phi;
}
double cxxGasPhase::Get_component_f(const std::string comp_name)
{
double f = -1.0;
for (size_t i = 0; i < this->gas_comps.size(); i++)
{
if (Utilities::strcmp_nocase(this->gas_comps[i].Get_phase_name().c_str(), comp_name.c_str()) == 0)
{
f = this->gas_comps[i].Get_f();
break;
}
}
return f;
}
cxxGasComp *
cxxGasPhase::Find_comp(const char * comp_name)
{

7
src/phreeqcpp/GasPhase.h
View File

@ -70,7 +70,12 @@ class cxxGasPhase:public cxxNumKeyword
cxxGasComp *Find_comp(const char * comp_name);
void Serialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles);
void Deserialize(Dictionary & dictionary, std::vector < int >&ints, std::vector < double >&doubles, int &ii, int &dd);
void Delete_component(const std::string comp_name);
void Set_component_moles(const std::string comp_name, const double moles);
double Get_component_moles(const std::string comp_name);
double Get_component_p(const std::string comp_name);
double Get_component_phi(const std::string comp_name);
double Get_component_f(const std::string comp_name);
protected:
void add(const cxxGasPhase & addee, LDBLE extensive);

80
src/phreeqcpp/Solution.cxx
View File

@ -83,6 +83,7 @@ cxxSolution::operator =(const cxxSolution &rhs)
this->isotopes = rhs.isotopes;
this->species_map = rhs.species_map;
this->log_gamma_map = rhs.log_gamma_map;
this->log_molalities_map = rhs.log_molalities_map;
if (this->initial_data)
delete initial_data;
if (rhs.initial_data != NULL)
@ -327,6 +328,19 @@ cxxSolution::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) con
s_oss << it->first << " " << it->second << "\n";
}
}
// log_molalities_map
if (log_molalities_map.size() > 0)
{
s_oss << indent1;
s_oss << "-log_molalities_map" << "\n";
std::map<int, double>::const_iterator it = this->log_molalities_map.begin();
for (; it != log_molalities_map.end(); it++)
{
s_oss << indent2;
s_oss << it->first << " " << it->second << "\n";
}
}
return;
}
@ -1013,7 +1027,32 @@ cxxSolution::read_raw(CParser & parser, bool check)
}
opt_save = CParser::OPT_DEFAULT;
break;
case 27: // log_molalities_map
{
int s_num;
if (parser.peek_token() != CParser::TT_EMPTY)
{
if (!(parser.get_iss() >> s_num))
{
parser.incr_input_error();
parser.error_msg("Expected integer for species number.",
PHRQ_io::OT_CONTINUE);
}
else
{
double d;
if (!(parser.get_iss() >> d))
{
parser.incr_input_error();
parser.error_msg("Expected double for species molality.",
PHRQ_io::OT_CONTINUE);
}
this->log_molalities_map[s_num] = d;
}
}
opt_save = 27;
}
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
@ -1415,6 +1454,19 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
this->log_gamma_map[git->first] = git->second;
}
}
// Add molalities
std::map<int, double>::const_iterator mit = addee.log_molalities_map.begin();
for (; mit != addee.log_molalities_map.end(); mit++)
{
if (this->log_molalities_map.find(mit->first) != this->log_molalities_map.end())
{
this->log_molalities_map[mit->first] = this->log_molalities_map[mit->first] * f1 + mit->second * f2;
}
else
{
this->log_molalities_map[mit->first] = mit->second;
}
}
}
}
@ -1618,6 +1670,18 @@ cxxSolution::Serialize(Dictionary & dictionary, std::vector < int >&ints,
doubles.push_back(it->second);
}
}
/*
* log_molalities_map
*/
ints.push_back((int)log_molalities_map.size());
{
std::map < int, double >::iterator it;
for (it = log_molalities_map.begin(); it != log_molalities_map.end(); it++)
{
ints.push_back(it->first);
doubles.push_back(it->second);
}
}
}
/* ---------------------------------------------------------------------- */
@ -1692,6 +1756,17 @@ cxxSolution::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
log_gamma_map[ints[ii++]] = doubles[dd++];
}
}
/*
* log_molalities_map
*/
{
log_molalities_map.clear();
int n = ints[ii++];
for (int i = 0; i < n; i++)
{
log_molalities_map[ints[ii++]] = doubles[dd++];
}
}
}
@ -1722,6 +1797,7 @@ const std::vector< std::string >::value_type temp_vopts[] = {
std::vector< std::string >::value_type("soln_vol"), // 23
std::vector< std::string >::value_type("species_map"), // 24
std::vector< std::string >::value_type("log_gamma_map"), // 25
std::vector< std::string >::value_type("potential") // 26
std::vector< std::string >::value_type("potential"), // 26
std::vector< std::string >::value_type("log_molalities_map") // 27
};
const std::vector< std::string > cxxSolution::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]);

2
src/phreeqcpp/Solution.h
View File

@ -66,6 +66,7 @@ class cxxSolution:public cxxNumKeyword
cxxNameDouble & Get_species_gamma(void) {return this->species_gamma;}
std::map<int, double> & Get_species_map(void) {return this->species_map;}
std::map<int, double> & Get_log_gamma_map(void) {return this->log_gamma_map;}
std::map<int, double>& Get_log_molalities_map(void) { return this->log_molalities_map; }
std::map < std::string, cxxSolutionIsotope > & Get_isotopes(void) {return this->isotopes;}
const std::map < std::string, cxxSolutionIsotope > & Get_isotopes(void)const {return this->isotopes;}
void Set_isotopes(const std::map < std::string, cxxSolutionIsotope > &iso ) {this->isotopes = iso;}
@ -141,6 +142,7 @@ class cxxSolution:public cxxNumKeyword
const static std::vector < std::string > vopts;
std::map<int, double> species_map;
std::map<int, double> log_gamma_map;
std::map<int, double> log_molalities_map;
};
#endif // !defined(SOLUTION_H_INCLUDED)

55
src/phreeqcpp/mainsubs.cpp
View File

@ -1352,19 +1352,19 @@ int Phreeqc::
xgas_save(int n_user)
/* ---------------------------------------------------------------------- */
{
/*
* Save gas composition into structure gas_phase with user
* number n_user.
*/
/*
* Save gas composition into structure gas_phase with user
* number n_user.
*/
char token[MAX_LENGTH];
if (use.Get_gas_phase_ptr() == NULL)
return (OK);
cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr();
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
cxxGasPhase temp_gas_phase(*gas_phase_ptr);
/*
* Store in gas_phase
*/
/*
* Store in gas_phase
*/
temp_gas_phase.Set_n_user(n_user);
temp_gas_phase.Set_n_user_end(n_user);
sprintf(token, "Gas phase after simulation %d.", simulation);
@ -1372,22 +1372,38 @@ xgas_save(int n_user)
temp_gas_phase.Set_new_def(false);
temp_gas_phase.Set_solution_equilibria(false);
temp_gas_phase.Set_n_solution(-99);
/*
* Update amounts
*/
for (size_t i = 0 ; i < temp_gas_phase.Get_gas_comps().size(); i++)
/*
* Update amounts
*/
bool PR = false;
if (gas_phase_ptr->Get_v_m() >= 0.01) PR = true;
for (size_t i = 0; i < temp_gas_phase.Get_gas_comps().size(); i++)
{
cxxGasComp * gc_ptr = &(temp_gas_phase.Get_gas_comps()[i]);
cxxGasComp* gc_ptr = &(temp_gas_phase.Get_gas_comps()[i]);
int k;
struct phase *phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE);
struct phase* phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE);
assert(phase_ptr);
gc_ptr->Set_moles(phase_ptr->moles_x);
if (PR)
{
gc_ptr->Set_moles(phase_ptr->moles_x);
gc_ptr->Set_p(phase_ptr->p_soln_x);
gc_ptr->Set_phi(phase_ptr->pr_phi);
gc_ptr->Set_f(phase_ptr->p_soln_x * phase_ptr->pr_phi);
}
else
{
gc_ptr->Set_moles(phase_ptr->moles_x);
gc_ptr->Set_p(phase_ptr->p_soln_x);
gc_ptr->Set_phi(1.0);
gc_ptr->Set_f(phase_ptr->p_soln_x);
}
}
Rxn_gas_phase_map[n_user] = temp_gas_phase;
use.Set_gas_phase_ptr(NULL);
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
xss_assemblage_save(int n_user)
@ -1671,6 +1687,15 @@ xsolution_save(int n_user)
temp_solution.Get_log_gamma_map()[s_x[i]->number] = s_x[i]->lg;
}
}
// saves molalities
temp_solution.Get_log_molalities_map().clear();
for (int i = 0; i < this->count_s_x; i++)
{
if (s_x[i]->type <= H2O)
{
temp_solution.Get_log_molalities_map()[s_x[i]->number] = s_x[i]->lm;
}
}
}
/*
* Save solution

36
src/phreeqcpp/tidy.cpp
View File

@ -971,9 +971,14 @@ tidy_gas_phase(void)
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
if (!PR)
gas_phase_ptr->Get_gas_comps()[j].Set_moles(
gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * gas_phase_ptr->Get_volume() /
R_LITER_ATM / gas_phase_ptr->Get_temperature());
{
double moles = gas_phase_ptr->Get_gas_comps()[j].Get_p_read() * gas_phase_ptr->Get_volume() /
R_LITER_ATM / gas_phase_ptr->Get_temperature();
gas_phase_ptr->Get_gas_comps()[j].Set_moles(moles);
gas_phase_ptr->Get_gas_comps()[j].Set_p(gas_phase_ptr->Get_gas_comps()[j].Get_p_read());
gas_phase_ptr->Get_gas_comps()[j].Set_phi(1.0);
gas_phase_ptr->Get_gas_comps()[j].Set_f(gas_phase_ptr->Get_gas_comps()[j].Get_p_read());
}
}
else
{
@ -999,10 +1004,15 @@ tidy_gas_phase(void)
{
P += gas_phase_ptr->Get_gas_comps()[j].Get_p_read();
if (!PR)
gas_phase_ptr->Get_gas_comps()[j].Set_moles (
gas_phase_ptr->Get_gas_comps()[j].Get_p_read() *
gas_phase_ptr->Get_volume() / R_LITER_ATM /
gas_phase_ptr->Get_temperature());
{
double moles = gas_phase_ptr->Get_gas_comps()[j].Get_p_read() *
gas_phase_ptr->Get_volume() / R_LITER_ATM /
gas_phase_ptr->Get_temperature();
gas_phase_ptr->Get_gas_comps()[j].Set_moles(moles);
gas_phase_ptr->Get_gas_comps()[j].Set_p(gas_phase_ptr->Get_gas_comps()[j].Get_p_read());
gas_phase_ptr->Get_gas_comps()[j].Set_phi(1.0);
gas_phase_ptr->Get_gas_comps()[j].Set_f(gas_phase_ptr->Get_gas_comps()[j].Get_p_read());
}
}
else
{
@ -1027,7 +1037,13 @@ tidy_gas_phase(void)
int k;
struct phase *phase_ptr = phase_bsearch(gas_phase_ptr->Get_gas_comps()[j_PR].Get_phase_name().c_str(), &k, FALSE);
if (gc[j_PR].Get_p_read() == 0)
{
gc[j_PR].Set_moles(0.0);
gc[j_PR].Set_p(0.0);
gc[j_PR].Set_phi(1.0);
gc[j_PR].Set_f(0.0);
continue;
}
if (phase_ptr)
{
phase_ptr->moles_x = gc[j_PR].Get_p_read() / P;
@ -1047,11 +1063,17 @@ tidy_gas_phase(void)
if (gc[j_PR].Get_p_read() == 0)
{
gc[j_PR].Set_moles(0.0);
gc[j_PR].Set_p(0.0);
gc[j_PR].Set_phi(1.0);
gc[j_PR].Set_f(0.0);
} else
{
if (phase_ptr)
{
gc[j_PR].Set_moles(phase_ptr->moles_x * gas_phase_ptr->Get_volume() / V_m);
gc[j_PR].Set_p(gc[j_PR].Get_p_read());
gc[j_PR].Set_phi(phase_ptr->pr_phi);
gc[j_PR].Set_f(gc[j_PR].Get_p_read()* phase_ptr->pr_phi);
gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + gc[j_PR].Get_moles());
}
}