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 '9d480bc124a4e9e979b4995476d086f90e3833e7'

Browse Source
This commit is contained in:
Darth Vader 2020-08-13 17:49:43 +00:00
commit fb6c36a792
5 changed files with 582 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
phreeqcpp/ExchComp.cxx
View File

@ -406,7 +406,7 @@ cxxExchComp::multiply(LDBLE extensive)
{
this->totals.multiply(extensive);
this->charge_balance *= extensive;
this->phase_proportion *= extensive;
//this->phase_proportion *= extensive;
}
void

2
phreeqcpp/PBasic.cpp
View File

@ -3993,7 +3993,7 @@ factor(struct LOC_exec * LINK)
// call callback Basic function
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->basic_callback(x1, x2, str);
PhreeqcPtr->PHRQ_free(str);
}
break;

4
phreeqcpp/Phreeqc.h
View File

@ -1033,15 +1033,19 @@ public:
int tidy_logk(void);
int tidy_exchange(void);
int tidy_min_exchange(void);
int update_min_exchange(void);
int tidy_kin_exchange(void);
int update_kin_exchange(void);
int tidy_gas_phase(void);
int tidy_inverse(void);
int tidy_isotopes(void);
int tidy_isotope_ratios(void);
int tidy_isotope_alphas(void);
int tidy_kin_surface(void);
int update_kin_surface(void);
int tidy_master_isotope(void);
int tidy_min_surface(void);
int update_min_surface(void);
int tidy_phases(void);
int tidy_pp_assemblage(void);
int tidy_solutions(void);

9
phreeqcpp/mainsubs.cpp
View File

@ -2679,6 +2679,11 @@ save_init(int i)
void
Phreeqc::do_mixes(void)
{
bool surf, exch, kin, min;
surf = (Rxn_surface_mix_map.size() > 0);
exch = (Rxn_exchange_mix_map.size() > 0);
kin = (Rxn_kinetics_mix_map.size() > 0);
min = (Rxn_pp_assemblage_mix_map.size() > 0);
Utilities::Rxn_mix(Rxn_solution_mix_map, Rxn_solution_map, this);
Utilities::Rxn_mix(Rxn_exchange_mix_map, Rxn_exchange_map, this);
Utilities::Rxn_mix(Rxn_gas_phase_mix_map, Rxn_gas_phase_map, this);
@ -2686,4 +2691,8 @@ Phreeqc::do_mixes(void)
Utilities::Rxn_mix(Rxn_pp_assemblage_mix_map, Rxn_pp_assemblage_map, this);
Utilities::Rxn_mix(Rxn_ss_assemblage_mix_map, Rxn_ss_assemblage_map, this);
Utilities::Rxn_mix(Rxn_surface_mix_map, Rxn_surface_map, this);
if (exch || kin) update_kin_exchange();
if (exch || min) update_min_exchange();
if (surf || min) update_min_surface();
if (surf || kin) update_kin_surface();
}

568
phreeqcpp/tidy.cpp
View File

@ -292,7 +292,45 @@ tidy_model(void)
{
tidy_solutions();
}
/*
* need to update exchange and surface related in case anything has changed
*/
if (keycount[Keywords::KEY_KINETICS] > 0 ||
keycount[Keywords::KEY_KINETICS_RAW] > 0 ||
keycount[Keywords::KEY_KINETICS_MODIFY] ||
keycount[Keywords::KEY_EXCHANGE] > 0 ||
keycount[Keywords::KEY_EXCHANGE_RAW] > 0 ||
keycount[Keywords::KEY_EXCHANGE_MODIFY])
{
update_kin_exchange();
}
if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] ||
keycount[Keywords::KEY_EXCHANGE] > 0 ||
keycount[Keywords::KEY_EXCHANGE_RAW] > 0 ||
keycount[Keywords::KEY_EXCHANGE_MODIFY])
{
update_min_exchange();
}
if (keycount[Keywords::KEY_EQUILIBRIUM_PHASES] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_RAW] > 0 ||
keycount[Keywords::KEY_EQUILIBRIUM_PHASES_MODIFY] ||
keycount[Keywords::KEY_SURFACE] > 0 ||
keycount[Keywords::KEY_SURFACE_RAW] > 0 ||
keycount[Keywords::KEY_SURFACE_MODIFY] > 0)
{
update_min_surface();
}
if (keycount[Keywords::KEY_KINETICS] > 0 ||
keycount[Keywords::KEY_KINETICS_RAW] > 0 ||
keycount[Keywords::KEY_KINETICS_MODIFY] > 0 ||
keycount[Keywords::KEY_SURFACE] > 0 ||
keycount[Keywords::KEY_SURFACE_RAW] > 0 ||
keycount[Keywords::KEY_SURFACE_MODIFY] > 0)
{
update_kin_surface();
}
/* if (new_model || new_exchange || new_pp_assemblage || new_surface || new_gas_phase || new_kinetics) reset_last_model(); */
if (new_model)
{
@ -3669,6 +3707,120 @@ tidy_kin_exchange(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_kin_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion. Exchange needs to be updated if the
* amount of kinetic reaction has changed. Corner case of
* zero moles.
*/
{
cxxKinetics* kinetics_ptr;
char* ptr;
LDBLE conc;
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
for ( ; it != Rxn_exchange_map.end(); it++)
{
cxxExchange* exchange_ptr = &(it->second);
if (exchange_ptr->Get_n_user() < 0) continue;
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
{
cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j];
if (comp_ref.Get_rate_name().size() == 0) continue;
double comp_moles = 0.0;
/* First find exchange master species */
cxxNameDouble nd = comp_ref.Get_totals();
cxxNameDouble::iterator kit = nd.begin();
bool found_exchange = false;
for (; kit != nd.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr == NULL || elt_ptr->master == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
kit->first.c_str());
error_msg(error_string, CONTINUE);
continue;
}
if (elt_ptr->master->type == EX)
{
comp_moles = kit->second;
found_exchange = true;
}
}
//if (!found_exchange)
//{
// input_error++;
// error_string = sformatf(
// "Exchange formula does not contain an exchange master species, %s",
// comp_ref.Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
/* Now find associated kinetic reaction ... */
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, exchange_ptr->Get_n_user())) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use exchange related to kinetic reaction, %s",
exchange_ptr->Get_n_user(), comp_ref.Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
size_t k;
for (k = 0; k < kinetics_ptr->Get_kinetics_comps().size(); k++)
{
if (strcmp_nocase
(comp_ref.Get_rate_name().c_str(),
kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to exchanger, %s, not found in KINETICS %d",
comp_ref.Get_rate_name().c_str(), comp_ref.Get_formula().c_str(), exchange_ptr->Get_n_user());
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
comp_ref.Set_rate_name(kinetics_ptr->Get_kinetics_comps()[k].Get_rate_name().c_str());
/* make exchanger concentration proportional to mineral ... */
conc = kinetics_ptr->Get_kinetics_comps()[k].Get_m() * comp_ref.Get_phase_proportion();
if (found_exchange && comp_moles > 0.0)
{
comp_ref.multiply(conc / comp_moles);
}
else /* need to generate totals from scratch */
{
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comp_ref.Set_totals(elt_list_NameDouble());
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_min_exchange(void)
/* ---------------------------------------------------------------------- */
/*
@ -3832,6 +3984,167 @@ tidy_min_exchange(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_min_exchange(void)
/* ---------------------------------------------------------------------- */
/*
* If exchanger is related to mineral, exchanger amount is
* set in proportion. Need to check in case exchange or min
* are modified.
*/
{
int n, jj;
char* ptr;
LDBLE conc;
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
for ( ; it != Rxn_exchange_map.end(); it++)
{
cxxExchange* exchange_ptr = &(it->second);
if (exchange_ptr->Get_n_user() < 0) continue;
n = exchange_ptr->Get_n_user();
// check elements
for (size_t j = 0; j < exchange_ptr->Get_exchange_comps().size(); j++)
{
double comp_moles = 0.0;
cxxExchComp& comp_ref = exchange_ptr->Get_exchange_comps()[j];
if (comp_ref.Get_phase_name().size() == 0) continue;
/* First find exchange master species */
cxxNameDouble nd = comp_ref.Get_totals();
cxxNameDouble::iterator kit = nd.begin();
bool found_exchange = false;
for (; kit != nd.end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
if (elt_ptr == NULL || elt_ptr->master == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
kit->first.c_str());
error_msg(error_string, CONTINUE);
continue;
}
if (elt_ptr->master->type == EX)
{
comp_moles = kit->second;
found_exchange = true;
}
}
//if (!found_exchange)
//{
// input_error++;
// error_string = sformatf(
// "Exchange formula does not contain an exchange master species, %s",
// comp_ref.Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n);
/* Now find the mineral on which exchanger depends... */
if (pp_assemblage_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Equilibrium_phases %d must be defined to use exchange related to mineral phase, %s",
n, comp_ref.Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
std::map<std::string, cxxPPassemblageComp>::iterator jit;
jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin();
for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++)
{
if (strcmp_nocase(comp_ref.Get_phase_name().c_str(), jit->first.c_str()) == 0)
{
break;
}
}
if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end())
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d",
comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
comp_ref.Set_phase_name(jit->first.c_str());
/* make exchanger concentration proportional to mineral ... */
conc = jit->second.Get_moles() * comp_ref.Get_phase_proportion();
if (found_exchange && comp_moles > 0.0)
{
comp_ref.multiply(conc / comp_moles);
}
else /* comp_moles is zero, need to redefine totals from scratch */
{
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
}
comp_ref.Set_totals(elt_list_NameDouble());
/*
* make sure exchange elements are in phase
*/
count_elts = 0;
paren_count = 0;
{
char* temp_formula = string_duplicate(comp_ref.Get_formula().c_str());
ptr = temp_formula;
get_elts_in_species(&ptr, -comp_ref.Get_phase_proportion());
free_check_null(temp_formula);
}
int l;
struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE);
if (phase_ptr != NULL)
{
char* temp_formula = string_duplicate(phase_ptr->formula);
ptr = temp_formula;
get_elts_in_species(&ptr, 1.0);
free_check_null(temp_formula);
}
else
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to exchanger, %s, not found in Equilibrium_Phases %d",
comp_ref.Get_phase_name().c_str(), comp_ref.Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
qsort(elt_list, (size_t)count_elts,
(size_t)sizeof(struct elt_list), elt_list_compare);
elt_list_combine();
for (jj = 0; jj < count_elts; jj++)
{
if (elt_list[jj].elt->primary->s->type != EX
&& elt_list[jj].coef < 0)
{
input_error++;
error_string = sformatf(
"Stoichiometry of exchanger, %s * %g mol sites/mol phase,\n\tmust be a subset of the related phase %s, %s.",
comp_ref.Get_formula().c_str(),
(double)comp_ref.Get_phase_proportion(),
phase_ptr->name,
phase_ptr->formula);
error_msg(error_string, CONTINUE);
break;
}
}
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_min_surface(void)
/* ---------------------------------------------------------------------- */
/*
@ -4084,6 +4397,134 @@ tidy_min_surface(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_min_surface(void)
/* ---------------------------------------------------------------------- */
/*
* If surface is related to mineral, surface amount is
* set in proportion
*/
{
std::map<int, cxxSurface>::iterator kit;
for (kit = Rxn_surface_map.begin(); kit != Rxn_surface_map.end(); kit++)
{
cxxSurface* surface_ptr = &(kit->second);
if (surface_ptr->Get_n_user() < 0) continue;
for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++)
{
double comp_moles = 0.0;
cxxSurfaceComp* surface_comp_ptr = &(surface_ptr->Get_surface_comps()[j]);
if (surface_comp_ptr->Get_phase_name().size() == 0) continue;
cxxSurfaceCharge* surface_charge_ptr = surface_ptr->Find_charge(surface_comp_ptr->Get_charge_name());
int n = surface_ptr->Get_n_user();
/* First find surface master species */
cxxNameDouble::iterator it;
for (it = surface_comp_ptr->Get_totals().begin(); it != surface_comp_ptr->Get_totals().end(); it++)
{
/* Find master species */
struct element* elt_ptr = element_store(it->first.c_str());
struct master* master_ptr = elt_ptr->master;
if (master_ptr == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
elt_ptr->name);
error_msg(error_string, CONTINUE);
continue;
}
if (master_ptr->type != SURF) continue;
comp_moles = it->second;
surface_comp_ptr->Set_master_element(elt_ptr->name);
break;
}
//if (surface_comp_ptr->Get_master_element().size() == 0)
//{
// input_error++;
// error_string = sformatf(
// "Surface formula does not contain a surface master species, %s",
// surface_comp_ptr->Get_formula().c_str());
// error_msg(error_string, CONTINUE);
// continue;
//}
/* Now find the mineral on which surface depends... */
cxxPPassemblage* pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, n);
if (pp_assemblage_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Equilibrium_phases %d must be defined to use surface related to mineral phase, %s",
n, surface_comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
std::map<std::string, cxxPPassemblageComp>::iterator jit;
jit = pp_assemblage_ptr->Get_pp_assemblage_comps().begin();
for (; jit != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); jit++)
{
if (strcmp_nocase(surface_comp_ptr->Get_phase_name().c_str(),
jit->first.c_str()) == 0)
{
break;
}
}
if (jit == pp_assemblage_ptr->Get_pp_assemblage_comps().end())
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to surface, %s, not found in Equilibrium_Phases %d",
surface_comp_ptr->Get_phase_name().c_str(), surface_comp_ptr->Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
int l;
struct phase* phase_ptr = phase_bsearch(jit->first.c_str(), &l, FALSE);
if (phase_ptr == NULL)
{
input_error++;
error_string = sformatf(
"Mineral, %s, related to surface, %s, not found in database.",
jit->first.c_str(), surface_comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
/* use database name for phase */
surface_comp_ptr->Set_phase_name(phase_ptr->name);
/* make surface concentration proportional to mineral ... */
LDBLE conc = jit->second.Get_moles() * surface_comp_ptr->Get_phase_proportion();
double grams = surface_charge_ptr->Get_grams();
if (comp_moles > 0.0)
{
surface_comp_ptr->multiply(conc / comp_moles);
}
else /* need to generate from scratch */
{
char* temp_formula = string_duplicate(surface_comp_ptr->Get_formula().c_str());
char* ptr = temp_formula;
count_elts = 0;
paren_count = 0;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
cxxNameDouble nd = elt_list_NameDouble();
surface_comp_ptr->Set_totals(nd);
}
if (grams > 0.0)
{
surface_charge_ptr->multiply(jit->second.Get_moles() / grams);
}
else /* need to generate from scratch */
{
surface_charge_ptr->Set_grams(jit->second.Get_moles());
surface_charge_ptr->Set_charge_balance(0.0);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
tidy_kin_surface(void)
/* ---------------------------------------------------------------------- */
/*
@ -4367,6 +4808,131 @@ tidy_kin_surface(void)
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
update_kin_surface(void)
/* ---------------------------------------------------------------------- */
/*
* If surface is related to mineral, surface amount is
* set in proportion. Need to update surface if
* moles of kinetic reaction changes
*/
{
cxxKinetics* kinetics_ptr;
std::map<int, cxxSurface>::iterator it;
for (it = Rxn_surface_map.begin(); it != Rxn_surface_map.end(); it++)
{
cxxSurface* surface_ptr = &(it->second);
if (surface_ptr->Get_n_user() < 0) continue;
int n = surface_ptr->Get_n_user();
for (size_t j = 0; j < surface_ptr->Get_surface_comps().size(); j++)
{
double comp_moles = 0.0;
cxxSurfaceComp* comp_ptr = &(surface_ptr->Get_surface_comps()[j]);
if (comp_ptr->Get_rate_name().size() == 0) continue;
comp_ptr->Set_master_element("");
/* First find surface master species */
int k;
cxxNameDouble::iterator kit;
for (kit = comp_ptr->Get_totals().begin(); kit != comp_ptr->Get_totals().end(); kit++)
{
/* Find master species */
struct element* elt_ptr = element_store(kit->first.c_str());
struct master* master_ptr = elt_ptr->master;
if (master_ptr == NULL)
{
input_error++;
error_string = sformatf("Master species not in database "
"for %s, skipping element.",
elt_ptr->name);
error_msg(error_string, CONTINUE);
continue;
}
if (master_ptr->type != SURF) continue;
comp_ptr->Set_master_element(elt_ptr->name);
comp_moles = kit->second;
break;
}
if (comp_ptr->Get_master_element().size() == 0)
{
input_error++;
error_string = sformatf(
"Surface formula does not contain a surface master species, %s",
comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
/* Now find the kinetic reaction on which surface depends... */
if ((kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, n)) == NULL)
{
input_error++;
error_string = sformatf(
"Kinetics %d must be defined to use surface related to kinetic reaction, %s",
n, comp_ptr->Get_formula().c_str());
error_msg(error_string, CONTINUE);
continue;
}
for (k = 0; k < (int)kinetics_ptr->Get_kinetics_comps().size(); k++)
{
cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
if (strcmp_nocase
(comp_ptr->Get_rate_name().c_str(),
kin_comp_ptr->Get_rate_name().c_str()) == 0)
{
break;
}
}
if (k == (int)kinetics_ptr->Get_kinetics_comps().size())
{
input_error++;
error_string = sformatf(
"Kinetic reaction, %s, related to surface, %s, not found in Kinetics %d",
comp_ptr->Get_rate_name().c_str(), comp_ptr->Get_formula().c_str(), n);
error_msg(error_string, CONTINUE);
continue;
}
cxxKineticsComp* kin_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[k]);
/* use database name for rate */
comp_ptr->Set_rate_name(kin_comp_ptr->Get_rate_name().c_str());
cxxSurfaceCharge* charge_ptr = surface_ptr->Find_charge(comp_ptr->Get_charge_name());
/* make surface concentration proportional to mineral ... */
LDBLE conc = kin_comp_ptr->Get_m() * comp_ptr->Get_phase_proportion();
double grams = charge_ptr->Get_grams();
if (comp_moles > 0.0)
{
comp_ptr->multiply(conc / comp_moles);
}
else /* need to generate from scratch */
{
char* temp_formula = string_duplicate(comp_ptr->Get_formula().c_str());
char* ptr = temp_formula;
count_elts = 0;
paren_count = 0;
get_elts_in_species(&ptr, conc);
free_check_null(temp_formula);
cxxNameDouble nd = elt_list_NameDouble();
comp_ptr->Set_totals(nd);
}
if (grams > 0.0)
{
charge_ptr->multiply(kin_comp_ptr->Get_m() / grams);
}
else /* need to generate from scratch */
{
charge_ptr->Set_grams(kin_comp_ptr->Get_m());
charge_ptr->Set_charge_balance(0.0);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
ss_prep(LDBLE t, cxxSS *ss_ptr, int print)
/* ---------------------------------------------------------------------- */
{