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

Working on ss

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/branches/ErrorHandling@6069 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2012-01-17 16:51:36 +00:00
parent 6bd23b9b8b
commit 2673f81bc7
10 changed files with 1501 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

431
SScomp.cxx Normal file
View File

@ -0,0 +1,431 @@
// SScomp.cxx: implementation of the cxxSScomp class.
//
//////////////////////////////////////////////////////////////////////
#ifdef _DEBUG
#pragma warning(disable : 4786) // disable truncation warning (Only used by debugger)
#endif
#include <cassert> // assert
#include <algorithm> // std::sort
#include "Utils.h" // define first
#include "Phreeqc.h"
#include "SScomp.h"
//#include "Dictionary.h"
#include "phqalloc.h"
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
cxxSScomp::cxxSScomp(PHRQ_io *io)
:
PHRQ_base(io)
//
// default constructor for cxxSScomp
//
{
name = "";
initial_moles = 0;
moles = 0;
init_moles = 0;
delta = 0;
fraction_x = 0;
log10_lambda = 0;
log10_fraction_x = 0;
dn = 0;
dnc = 0;
dnb = 0;
}
#ifdef SKIP
cxxSScomp::cxxSScomp(struct pure_phase * pure_phase_ptr, PHRQ_io *io)
:
PHRQ_base(io)
//
// constructor for cxxSScomp from struct pure_phase
//
{
this->Set_name(pure_phase_ptr->name);
this->Set_add_formula(pure_phase_ptr->add_formula);
si = pure_phase_ptr->si;
si_org = pure_phase_ptr->si_org;
moles = pure_phase_ptr->moles;
delta = pure_phase_ptr->delta;
initial_moles = pure_phase_ptr->initial_moles;
force_equality = (pure_phase_ptr->force_equality == TRUE);
dissolve_only = (pure_phase_ptr->dissolve_only == TRUE);
precipitate_only = (pure_phase_ptr->precipitate_only == TRUE);
}
#endif
cxxSScomp::~cxxSScomp()
{
}
void
cxxSScomp::dump_raw(std::ostream & s_oss, unsigned int indent) const
{
//const char ERR_MESSAGE[] = "Packing pure_phase message: %s, element not found\n";
unsigned int i;
s_oss.precision(DBL_DIG - 1);
std::string indent0(""), indent1(""), indent2("");
for (i = 0; i < indent; ++i)
indent0.append(Utilities::INDENT);
for (i = 0; i < indent + 1; ++i)
indent1.append(Utilities::INDENT);
for (i = 0; i < indent + 2; ++i)
indent2.append(Utilities::INDENT);
// Pure_Phase element and attributes
if (this->name.size() != 0)
s_oss << indent0 << "-name " << this->name << "\n";
//if (this->add_formula.size() != 0)
// s_oss << indent1 << "-add_formula " << this->add_formula << "\n";
s_oss << indent1 << "-initial_moles " << this->initial_moles << "\n";
s_oss << indent1 << "-moles " << this->moles << "\n";
s_oss << indent1 << "-initial_moles " << this->init_moles << "\n";
s_oss << indent1 << "-delta " << this->delta << "\n";
s_oss << indent1 << "-fraction_x " << this->fraction_x << "\n";
s_oss << indent1 << "-log10_lambda " << this->log10_lambda << "\n";
s_oss << indent1 << "-log10_fraction_x " << this->log10_fraction_x << "\n";
s_oss << indent1 << "-dn " << this->dn << "\n";
s_oss << indent1 << "-dnc " << this->dnc << "\n";
s_oss << indent1 << "-dnb " << this->dnb << "\n";
}
void
cxxSScomp::read_raw(CParser & parser, bool check)
{
std::string str;
static std::vector < std::string > vopts;
if (vopts.empty())
{
vopts.reserve(10);
vopts.push_back("name"); // 0
vopts.push_back("initial_moles"); // 1
vopts.push_back("moles"); // 2
vopts.push_back("init_moles"); // 3
vopts.push_back("delta"); // 4
vopts.push_back("fraction_x"); // 5
vopts.push_back("log10_lambda"); // 6
vopts.push_back("log10_fraction_x"); // 7
vopts.push_back("dn"); // 8
vopts.push_back("dnc"); // 9
vopts.push_back("dnb"); // 9
}
std::istream::pos_type ptr;
std::istream::pos_type next_char;
std::string token;
int opt_save;
opt_save = CParser::OPT_ERROR;
bool name_defined(false);
bool initial_moles_defined(false);
bool moles_defined(false);
//bool init_moles_defined(false);
//bool delta_defined(false);
for (;;)
{
int opt = parser.get_option(vopts, next_char);
if (opt == CParser::OPT_DEFAULT)
{
opt = opt_save;
}
switch (opt)
{
case CParser::OPT_EOF:
break;
case CParser::OPT_KEYWORD:
break;
case CParser::OPT_DEFAULT:
case CParser::OPT_ERROR:
opt = CParser::OPT_KEYWORD;
// Allow return to Exchange for more processing
//parser.error_msg("Unknown input in PURE_PHASE read.", PHRQ_io::OT_CONTINUE);
//parser.error_msg(parser.line().c_str(), PHRQ_io::OT_CONTINUE);
break;
case 0: // name
if (!(parser.get_iss() >> str))
{
this->name.clear();
parser.incr_input_error();
parser.error_msg("Expected string value for name.",
PHRQ_io::OT_CONTINUE);
}
else
{
this->name = str;
}
name_defined = true;
break;
case 1: // initial_moles
if (!(parser.get_iss() >> this->initial_moles))
{
this->initial_moles = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for initial_moles.",
PHRQ_io::OT_CONTINUE);
}
initial_moles_defined = true;
break;
case 2: // moles
if (!(parser.get_iss() >> this->moles))
{
this->moles = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for moles.",
PHRQ_io::OT_CONTINUE);
}
moles_defined = true;
break;
case 3: // init_moles
if (!(parser.get_iss() >> this->init_moles))
{
this->init_moles = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for init_moles.",
PHRQ_io::OT_CONTINUE);
}
//init_moles_defined = true;
break;
case 4: // delta
if (!(parser.get_iss() >> this->delta))
{
this->delta = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for delta.",
PHRQ_io::OT_CONTINUE);
}
//delta_defined = true;
break;
case 5: // fraction_x
if (!(parser.get_iss() >> this->fraction_x))
{
this->fraction_x = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for fraction_x.",
PHRQ_io::OT_CONTINUE);
}
//initial_moles_defined = true;
break;
case 6: // log10_lambda
if (!(parser.get_iss() >> this->log10_lambda))
{
this->log10_lambda = 0;
parser.incr_input_error();
parser.error_msg("Expected boolean value for log10_lambda.",
PHRQ_io::OT_CONTINUE);
}
break;
case 7: // log10_fraction_x
if (!(parser.get_iss() >> this->log10_fraction_x))
{
this->log10_fraction_x = 0;
parser.incr_input_error();
parser.error_msg("Expected boolean value for log10_fraction_x.",
PHRQ_io::OT_CONTINUE);
}
break;
case 8: // dn
if (!(parser.get_iss() >> this->dn))
{
this->dn = 0;
parser.incr_input_error();
parser.error_msg("Expected boolean value for dn.",
PHRQ_io::OT_CONTINUE);
}
break;
case 9: // dnc
if (!(parser.get_iss() >> this->dnc))
{
this->dnc = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for dnc.",
PHRQ_io::OT_CONTINUE);
}
//dnc_defined = true;
break;
case 10: // dnb
if (!(parser.get_iss() >> this->dnb))
{
this->dnb = 0;
parser.incr_input_error();
parser.error_msg("Expected numeric value for dnb.",
PHRQ_io::OT_CONTINUE);
}
//dnc_defined = true;
break;
}
if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD)
break;
}
// members that must be defined
if (check)
{
if (name_defined == false)
{
parser.incr_input_error();
parser.error_msg("Name not defined for PPassemblageComp input.",
PHRQ_io::OT_CONTINUE);
}
if (moles_defined == false)
{
parser.incr_input_error();
parser.error_msg("Moles not defined for PPassemblageComp input.",
PHRQ_io::OT_CONTINUE);
}
if (initial_moles_defined == false)
{
parser.incr_input_error();
parser.error_msg("Initial_moles not defined for PPassemblageComp input.",
PHRQ_io::OT_CONTINUE);
}
}
}
#ifdef USE_MPI
void
cxxSScomp::mpi_pack(std::vector < int >&ints,
std::vector < double >&doubles)
{
extern cxxDictionary dictionary;
ints.push_back(dictionary.string2int(this->name));
ints.push_back(dictionary.string2int(this->add_formula));
doubles.push_back(this->si);
doubles.push_back(this->si_org);
doubles.push_back(this->moles);
doubles.push_back(this->delta);
doubles.push_back(this->initial_moles);
ints.push_back((int) this->force_equality);
ints.push_back((int) this->dissolve_only);
ints.push_back((int) this->precipitate_only);
}
void
cxxSScomp::mpi_unpack(int *ints, int *ii, double *doubles, int *dd)
{
extern cxxDictionary dictionary;
int i = *ii;
int d = *dd;
this->name = dictionary.int2stdstring(ints[i++]);
this->add_formula = dictionary.int2stdstring(ints[i++]);
this->si = doubles[d++];
this->si_org = doubles[d++];
this->moles = doubles[d++];
this->delta = doubles[d++];
this->initial_moles = doubles[d++];
this->force_equality = (ints[i++] != 0);
this->dissolve_only = (ints[i++] != 0);
this->precipitate_only = (ints[i++] != 0);
*ii = i;
*dd = d;
}
#endif
#ifdef SKIP
void
cxxSScomp::totalize(PHREEQC_PTR_ARG)
{
this->totals.clear();
// component structures
if (this->add_formula.size() != 0)
return;
struct phase *phase_ptr;
int l;
phase_ptr = P_INSTANCE_POINTER phase_bsearch(this->name.c_str(), &l, FALSE);
if (phase_ptr != NULL)
{
cxxNameDouble phase_formula(phase_ptr->next_elt);
this->totals.add_extensive(phase_formula, this->moles);
}
else
{
assert(false);
}
return;
}
#endif
#ifdef SKIP
void
cxxSScomp::add(const cxxSScomp & addee, double extensive)
{
double ext1, ext2, f1, f2;
if (extensive == 0.0)
return;
if (addee.name.size() == 0)
return;
// this and addee must have same name
// otherwise generate a new PPassemblagComp with multiply
ext1 = this->moles;
ext2 = addee.moles * extensive;
if (ext1 + ext2 != 0)
{
f1 = ext1 / (ext1 + ext2);
f2 = ext2 / (ext1 + ext2);
}
else
{
f1 = 0.5;
f2 = 0.5;
}
//char * name;
//char *add_formula;
if (this->name.size() == 0 && addee.name.size() == 0)
{
return;
}
assert(this->name == addee.name);
if (this->add_formula != addee.add_formula)
{
std::ostringstream oss;
oss <<
"Can not mix two Equilibrium_phases with differing add_formulae., "
<< this->name;
error_msg(oss.str().c_str(), CONTINUE);
return;
}
//double si;
this->si = this->si * f1 + addee.si * f2;
//double si_org;
this->si_org = this->si_org * f1 + addee.si_org * f2;
//double moles;
this->moles += addee.moles * extensive;
//double delta;
this->delta += addee.delta * extensive;
//double initial_moles;
this->initial_moles += addee.initial_moles * extensive;
//bool force_equality;
//bool dissolve_only;
}
#endif
void
cxxSScomp::multiply(double extensive)
{
//char * name;
//char *add_formula;
//double si;
//double moles;
this->moles *= extensive;
//double delta;
this->delta *= extensive;
//double initial_moles;
this->initial_moles *= extensive;
//bool force_equality;
//bool dissolve_only;
}

67
SScomp.h Normal file
View File

@ -0,0 +1,67 @@
#if !defined(SSCOMP_H_INCLUDED)
#define SSCOMP_H_INCLUDED
#include <cassert> // assert
#include <map> // std::map
#include <string> // std::string
#include <list> // std::list
#include <vector> // std::vector
#include "NameDouble.h"
#include "Phreeqc_class.h"
class cxxSScomp: public PHRQ_base
{
public:
cxxSScomp(PHRQ_io *io=NULL);
~cxxSScomp();
void dump_raw(std::ostream & s_oss, unsigned int indent) const;
void read_raw(CParser & parser, bool check = true);
const std::string &Get_name() const {return this->name;}
void Set_name(std::string s) {this->name = s;}
LDBLE Get_initial_moles() const {return this->initial_moles;}
void Set_initial_moles(LDBLE t) {this->initial_moles = t;}
LDBLE Get_moles() const {return this->moles;}
void Set_moles(LDBLE t) {this->moles = t;}
LDBLE Get_init_moles() const {return this->init_moles;}
void Set_init_moles(LDBLE t) {this->init_moles = t;}
LDBLE Get_delta() const {return this->delta;}
void Set_delta(LDBLE t) {this->delta = t;}
LDBLE Get_fraction_x() const {return this->fraction_x;}
void Set_fraction_x(LDBLE t) {this->fraction_x = t;}
LDBLE Get_log10_lambda() const {return this->log10_lambda;}
void Set_log10_lambda(LDBLE t) {this->log10_lambda = t;}
LDBLE Get_log10_fraction_x() const {return this->log10_fraction_x;}
void Set_log10_fraction_x(LDBLE t) {this->log10_fraction_x = t;}
LDBLE Get_dn() const {return this->dn;}
void Set_dn(LDBLE t) {this->dn = t;}
LDBLE Get_dnc() const {return this->dnc;}
void Set_dnc(LDBLE t) {this->dnc = t;}
LDBLE Get_dnb() const {return this->dnb;}
void Set_dnb(LDBLE t) {this->dnb = t;}
//void add(const cxxSScomp & comp, double extensive);
void multiply(double extensive);
#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
protected:
std::string name;
//struct phase *phase;
LDBLE initial_moles;
LDBLE moles;
LDBLE init_moles;
LDBLE delta;
LDBLE fraction_x;
LDBLE log10_lambda;
LDBLE log10_fraction_x;
LDBLE dn, dnc, dnb;
};
#endif // !defined(SSCOMP_H_INCLUDED)

2
phreeqc/basicsubs.cpp
View File

@ -2794,7 +2794,7 @@ system_total_elt_secondary(const char *total_name)
/*
* Look for element
*/
for (int j1 = 0; j1 < (size_t) count_elts; j1++)
for (size_t j1 = 0; j1 < (size_t) count_elts; j1++)
{
if (strcmp(elt_list[j1].elt->name, total_name) == 0)
{

14
phreeqc/global_structures.h
View File

@ -84,7 +84,7 @@
#define SURFACE_CB1 22
#define SURFACE_CB2 23
#define GAS_MOLES 24
#define S_S_MOLES 25
#define SS_MOLES 25
#define PITZER_GAMMA 26
#define SLACK 28
/* state */
@ -648,6 +648,7 @@ struct name_coef
const char *name;
LDBLE coef;
};
#ifdef SKIP
/*----------------------------------------------------------------------
* Solid solution
*---------------------------------------------------------------------- */
@ -692,6 +693,7 @@ struct s_s_comp
LDBLE log10_fraction_x;
LDBLE dn, dnc, dnb;
};
#endif
#ifdef SKIP
/*----------------------------------------------------------------------
* Pure-phase assemblage
@ -1073,10 +1075,12 @@ struct unknown
const char * exch_comp;
//struct pure_phase *pure_phase;
const char *pp_assemblage_comp_name;
struct s_s *s_s;
struct s_s_comp *s_s_comp;
int s_s_comp_number;
int s_s_in;
//struct s_s *s_s;
//struct s_s_comp *s_s_comp;
const char * ss_name;
const char * ss_comp_name;
int ss_comp_number;
int ss_in;
struct surface_comp *surface_comp;
LDBLE related_moles;
struct unknown *potential_unknown, *potential_unknown1,

28
phreeqc/model.cpp
View File

@ -579,7 +579,7 @@ check_residuals(void)
x[i]->description, (double) residual[i]);
}
}
else if (x[i]->type == S_S_MOLES)
else if (x[i]->type == SS_MOLES)
{
if (x[i]->s_s_in == FALSE)
continue;
@ -1388,7 +1388,7 @@ ineq(int in_kode)
* Solid solution
*/
}
else if (x[i]->type == S_S_MOLES && x[i]->s_s_in == TRUE)
else if (x[i]->type == SS_MOLES && x[i]->s_s_in == TRUE)
{
memcpy((void *) &(ineq_array[l_count_rows * max_column_count]),
(void *) &(array[i * (count_unknowns + 1)]),
@ -1426,7 +1426,7 @@ ineq(int in_kode)
}
if (x[i]->type != SOLUTION_PHASE_BOUNDARY &&
x[i]->type != ALK &&
x[i]->type != GAS_MOLES && x[i]->type != S_S_MOLES
x[i]->type != GAS_MOLES && x[i]->type != SS_MOLES
/* && x[i]->type != PP */
&& x[i]->type != SLACK
)
@ -1773,7 +1773,7 @@ ineq(int in_kode)
{
for (i = s_s_unknown->number; i < count_unknowns; i++)
{
if (x[i]->type != S_S_MOLES)
if (x[i]->type != SS_MOLES)
break;
if (x[i]->phase->in == TRUE && x[i]->s_s_in == TRUE)
{
@ -1835,7 +1835,7 @@ ineq(int in_kode)
*/
for (i = 0; i < count_unknowns; i++)
{
if ((x[i]->type == PP || x[i]->type == S_S_MOLES)
if ((x[i]->type == PP || x[i]->type == SS_MOLES)
&& pp_column_scale != 1.0)
{
for (j = 0; j < l_count_rows; j++)
@ -2415,7 +2415,7 @@ mb_s_s(void)
}
for (i = s_s_unknown->number; i < count_unknowns; i++)
{
if (x[i]->type != S_S_MOLES)
if (x[i]->type != SS_MOLES)
break;
x[i]->s_s_in = x[i]->s_s->s_s_in;
}
@ -3108,7 +3108,7 @@ reset(void)
*/
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type == PP || x[i]->type == S_S_MOLES)
if (x[i]->type == PP || x[i]->type == SS_MOLES)
{
if (delta[i] < -1e8)
@ -3249,7 +3249,7 @@ reset(void)
for (i = 0; i < count_unknowns; i++)
{
x[i]->delta /= factor;
if (x[i]->type == PP || x[i]->type == S_S_MOLES)
if (x[i]->type == PP || x[i]->type == SS_MOLES)
delta[i] /= factor;
}
@ -3322,7 +3322,7 @@ reset(void)
up = 1.;
down = x[i]->moles;
}
else if (x[i]->type == S_S_MOLES)
else if (x[i]->type == SS_MOLES)
{
continue;
}
@ -3397,7 +3397,7 @@ reset(void)
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type != PP && x[i]->type != S_S_MOLES)
if (x[i]->type != PP && x[i]->type != SS_MOLES)
delta[i] *= factor;
}
@ -3774,7 +3774,7 @@ reset(void)
}
}
else if (x[i]->type == S_S_MOLES)
else if (x[i]->type == SS_MOLES)
{
/*if (fabs(delta[i]) > epsilon) converge=FALSE; */
@ -4115,7 +4115,7 @@ residuals(void)
converge = FALSE;
}
}
else if (x[i]->type == S_S_MOLES)
else if (x[i]->type == SS_MOLES)
{
residual[i] = x[i]->f * LOG_10;
//if (x[i]->moles <= MIN_TOTAL_SS && iterations > 2)
@ -5361,7 +5361,7 @@ numerical_jacobian(void)
reset();
d2 = delta[i];
break;
case S_S_MOLES:
case SS_MOLES:
if (x[i]->s_s_in == FALSE)
continue;
for (j = 0; j < count_unknowns; j++)
@ -5442,7 +5442,7 @@ numerical_jacobian(void)
delta[i] = -d2;
reset();
break;
case S_S_MOLES:
case SS_MOLES:
delta[i] = -d2;
reset();
break;

2
phreeqc/pitzer.cpp
View File

@ -1733,7 +1733,7 @@ Restart:
break;
case MU:
case PP:
case S_S_MOLES:
case SS_MOLES:
continue;
break;
}

379
phreeqc/prep.cpp
View File

@ -6,6 +6,8 @@
#include "Exchange.h"
#include "GasPhase.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "SS.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
prep(void)
@ -260,7 +262,40 @@ quick_setup(void)
{
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type == S_S_MOLES)
if (x[i]->type == SS_MOLES)
break;
}
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t j = 0; j < ss_ptrs.size(); j++)
{
for (size_t k = 0; k < ss_ptrs[j]->Get_ss_comps().size(); k++)
{
cxxSScomp *comp_ptr = &(ss_ptrs[j]->Get_ss_comps()[k]);
x[i]->moles = comp_ptr->Get_moles();
if (x[i]->moles <= 0)
{
x[i]->moles = MIN_TOTAL_SS;
comp_ptr->Set_moles(MIN_TOTAL_SS);
}
comp_ptr->Set_initial_moles(x[i]->moles);
x[i]->ln_moles = log(x[i]->moles);
x[i]->phase->dn = comp_ptr->Get_dn();
x[i]->phase->dnb = comp_ptr->Get_dnb();
x[i]->phase->dnc = comp_ptr->Get_dnc();
x[i]->phase->log10_fraction_x = comp_ptr->Get_log10_fraction_x();
x[i]->phase->log10_lambda = comp_ptr->Get_log10_lambda();
i++;
}
}
}
#ifdef SKIP
if (s_s_unknown != NULL)
{
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type == SS_MOLES)
break;
}
for (j = 0; j < use.Get_ss_assemblage_ptr()->count_s_s; j++)
@ -289,6 +324,7 @@ quick_setup(void)
}
}
}
#endif
/*
* exchange
*/
@ -666,6 +702,223 @@ int Phreeqc::
build_ss_assemblage(void)
/* ---------------------------------------------------------------------- */
{
/*
* Put coefficients into lists to sum iaps to test for equilibrium
* Put coefficients into lists to build jacobian for
* mass action equation for component
* mass balance equations for elements contained in solid solutions
*/
//int i, j, k, l, stop;
bool stop;
int row, col;
struct master *master_ptr;
struct rxn_token *rxn_ptr;
//struct s_s *s_s_ptr, *s_s_ptr_old;
//char token[MAX_LENGTH];
char *ptr;
if (s_s_unknown == NULL)
return (OK);
cxxSS * ss_ptr_old = NULL;
col = 0;
for (int i = 0; i < count_unknowns; i++)
{
if (x[i]->type != SS_MOLES)
continue;
cxxSS *ss_ptr = use.Get_ss_assemblage_ptr()->Find(x[i]->ss_name);
assert(ss_ptr);
if (ss_ptr != ss_ptr_old)
{
col = x[i]->number;
ss_ptr_old = ss_ptr;
}
/*
* Calculate function value (inverse saturation index)
*/
if (x[i]->phase->rxn_x == NULL)
continue;
store_mb(&(x[i]->phase->lk), &(x[i]->f), 1.0);
for (rxn_ptr = x[i]->phase->rxn_x->token + 1; rxn_ptr->s != NULL;
rxn_ptr++)
{
store_mb(&(rxn_ptr->s->la), &(x[i]->f), -rxn_ptr->coef);
}
/* include mole fraction */
store_mb(&(x[i]->phase->log10_fraction_x), &(x[i]->f), 1.0);
/* include activity coeficient */
store_mb(&(x[i]->phase->log10_lambda), &(x[i]->f), 1.0);
/*
* Put coefficients into mass action equations
*/
/* first IAP terms */
for (rxn_ptr = x[i]->phase->rxn_x->token + 1; rxn_ptr->s != NULL;
rxn_ptr++)
{
if (rxn_ptr->s->secondary != NULL
&& rxn_ptr->s->secondary->in == TRUE)
{
master_ptr = rxn_ptr->s->secondary;
}
else
{
master_ptr = rxn_ptr->s->primary;
}
if (master_ptr == NULL || master_ptr->unknown == NULL)
continue;
store_jacob0(x[i]->number, master_ptr->unknown->number,
rxn_ptr->coef);
}
if (ss_ptr->Get_a0() != 0.0 || ss_ptr->Get_a1() != 0.0)
{
/*
* For binary solid solution
*/
/* next dnc terms */
row = x[i]->number * (count_unknowns + 1);
if (x[i]->ss_comp_number == 0)
{
col = x[i]->number;
}
else
{
col = x[i]->number - 1;
}
store_jacob(&(x[i]->phase->dnc), &(array[row + col]), -1);
/* next dnb terms */
col++;
store_jacob(&(x[i]->phase->dnb), &(array[row + col]), -1);
}
else
{
/*
* For ideal solid solution
*/
row = x[i]->number * (count_unknowns + 1);
for (size_t j = 0; j < ss_ptr->Get_ss_comps().size(); j++)
{
if (j != x[i]->ss_comp_number)
{
/* store_jacob (&(s_s_ptr->dn), &(array[row + col + j]), -1.0); */
store_jacob(&(x[i]->phase->dn), &(array[row + col + j]),
-1.0);
}
else
{
store_jacob(&(x[i]->phase->dnb), &(array[row + col + j]),
-1.0);
}
}
}
/*
* Put coefficients into mass balance equations
*/
count_elts = 0;
paren_count = 0;
char * token = string_duplicate(x[i]->phase->formula);
ptr = token;
get_elts_in_species(&ptr, 1.0);
free_check_null(token);
/*
* Go through elements in phase
*/
#ifdef COMBINE
change_hydrogen_in_elt_list(0);
#endif
for (int j = 0; j < count_elts; j++)
{
if (strcmp(elt_list[j].elt->name, "H") == 0
&& mass_hydrogen_unknown != NULL)
{
store_jacob0(mass_hydrogen_unknown->number, x[i]->number,
-elt_list[j].coef);
store_sum_deltas(&(delta[i]), &mass_hydrogen_unknown->delta,
elt_list[j].coef);
}
else if (strcmp(elt_list[j].elt->name, "O") == 0
&& mass_oxygen_unknown != NULL)
{
store_jacob0(mass_oxygen_unknown->number, x[i]->number,
-elt_list[j].coef);
store_sum_deltas(&(delta[i]), &mass_oxygen_unknown->delta,
elt_list[j].coef);
}
else
{
master_ptr = elt_list[j].elt->primary;
if (master_ptr->in == FALSE)
{
master_ptr = master_ptr->s->secondary;
}
if (master_ptr == NULL || master_ptr->in == FALSE)
{
if (state != ADVECTION && state != TRANSPORT
&& state != PHAST)
{
error_string = sformatf(
"Element in phase, %s, is not in model.",
x[i]->phase->name);
warning_msg(error_string);
}
if (master_ptr != NULL)
{
master_ptr->s->la = -999.9;
}
/*
* Master species is in model
*/
}
else if (master_ptr->in == TRUE)
{
store_jacob0(master_ptr->unknown->number, x[i]->number,
-elt_list[j].coef);
store_sum_deltas(&delta[i], &master_ptr->unknown->delta,
elt_list[j].coef);
/*
* Master species in equation needs to be rewritten
*/
}
else if (master_ptr->in == REWRITE)
{
stop = FALSE;
for (int k = 0; k < count_unknowns; k++)
{
if (x[k]->type != MB)
continue;
for (int l = 0; x[k]->master[l] != NULL; l++)
{
if (x[k]->master[l] == master_ptr)
{
store_jacob0(x[k]->master[0]->unknown->
number, x[i]->number,
-elt_list[j].coef);
store_sum_deltas(&delta[i],
&x[k]->master[0]->unknown->
delta, elt_list[j].coef);
stop = TRUE;
break;
}
}
if (stop == TRUE)
break;
}
}
}
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
build_ss_assemblage(void)
/* ---------------------------------------------------------------------- */
{
/*
* Put coefficients into lists to sum iaps to test for equilibrium
* Put coefficients into lists to build jacobian for
@ -686,7 +939,7 @@ build_ss_assemblage(void)
col = 0;
for (i = 0; i < count_unknowns; i++)
{
if (x[i]->type != S_S_MOLES)
if (x[i]->type != SS_MOLES)
continue;
s_s_ptr = x[i]->s_s;
if (s_s_ptr != s_s_ptr_old)
@ -874,7 +1127,7 @@ build_ss_assemblage(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
build_jacobian_sums(int k)
@ -3488,6 +3741,58 @@ int Phreeqc::
setup_ss_assemblage(void)
/* ---------------------------------------------------------------------- */
{
/*
* Fill in data for solid solution unknowns (sum of partial pressures)
* in unknown structure
*/
//int i, j;
if (use.Get_ss_assemblage_ptr() == NULL)
return (OK);
/*
* One for each component in each solid solution
*/
s_s_unknown = NULL;
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t j = 0; j < ss_ptrs.size(); j++)
{
for (size_t i = 0; i < ss_ptrs[j]->Get_ss_comps().size(); i++)
{
cxxSScomp *comp_ptr = &(ss_ptrs[j]->Get_ss_comps()[i]);
int l;
struct phase* phase_ptr = phase_bsearch(comp_ptr->Get_name().c_str(), &l, FALSE);
x[count_unknowns]->type = SS_MOLES;
x[count_unknowns]->description = string_hsave(comp_ptr->Get_name().c_str());
x[count_unknowns]->moles = 0.0;
if (comp_ptr->Get_moles() <= 0)
{
comp_ptr->Set_moles(MIN_TOTAL_SS);
}
x[count_unknowns]->moles = comp_ptr->Get_moles();
comp_ptr->Set_initial_moles(x[count_unknowns]->moles);
x[count_unknowns]->ln_moles = log(x[count_unknowns]->moles);
x[count_unknowns]->ss_name = string_hsave(ss_ptrs[j]->Get_name().c_str());
x[count_unknowns]->ss_comp_name = string_hsave(comp_ptr->Get_name().c_str());
x[count_unknowns]->ss_comp_number = (int) i;
x[count_unknowns]->phase = phase_ptr;
x[count_unknowns]->number = count_unknowns;
x[count_unknowns]->phase->dn = comp_ptr->Get_dn();
x[count_unknowns]->phase->dnb = comp_ptr->Get_dnb();
x[count_unknowns]->phase->dnc = comp_ptr->Get_dnc();
x[count_unknowns]->phase->log10_fraction_x = comp_ptr->Get_log10_fraction_x();
x[count_unknowns]->phase->log10_lambda =comp_ptr->Get_log10_lambda();
if (s_s_unknown == NULL)
s_s_unknown = x[count_unknowns];
count_unknowns++;
}
}
return (OK);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
setup_ss_assemblage(void)
/* ---------------------------------------------------------------------- */
{
/*
* Fill in data for solid solution unknowns (sum of partial pressures)
* in unknown structure
@ -3503,7 +3808,7 @@ setup_ss_assemblage(void)
{
for (i = 0; i < use.Get_ss_assemblage_ptr()->s_s[j].count_comps; i++)
{
x[count_unknowns]->type = S_S_MOLES;
x[count_unknowns]->type = SS_MOLES;
x[count_unknowns]->description =
string_hsave(use.Get_ss_assemblage_ptr()->s_s[j].comps[i].name);
x[count_unknowns]->moles = 0.0;
@ -3540,7 +3845,7 @@ setup_ss_assemblage(void)
}
return (OK);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
setup_surface(void)
@ -4926,6 +5231,15 @@ setup_unknowns(void)
/*
* Count solid solutions
*/
if (use.Get_ss_assemblage_ptr() != NULL)
{
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t i = 0; i < ss_ptrs.size(); i++)
{
max_unknowns += ss_ptrs[i]->Get_ss_comps().size();
}
}
#ifdef SKIP
if (use.Get_ss_assemblage_ptr() != NULL)
{
/* max_unknowns += 2 * use.Get_ss_assemblage_ptr()->count_s_s; */
@ -4934,6 +5248,7 @@ setup_unknowns(void)
max_unknowns += use.Get_ss_assemblage_ptr()->s_s[i].count_comps;
}
}
#endif
/*
* Pitzer/Sit
@ -5751,6 +6066,18 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
/*
* Calculate miscibility gaps for solid solutions
*/
if (use.Get_ss_assemblage_ptr() != NULL)
{
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t i = 0; i < ss_ptrs.size(); i++)
{
if (fabs(tempk - ss_ptrs[i]->Get_tk()) > 0.01)
{
ss_prep(tempk, ss_ptrs[i], FALSE);
}
}
}
#ifdef SKIP
if (use.Get_ss_assemblage_ptr() != NULL)
{
for (i = 0; i < use.Get_ss_assemblage_ptr()->count_s_s; i++)
@ -5761,6 +6088,7 @@ k_temp(LDBLE tc, LDBLE pa) /* pa - pressure in atm */
}
}
}
#endif
return (OK);
}
@ -5865,6 +6193,26 @@ save_model(void)
last_model.ss_assemblage =
(const char **) free_check_null(last_model.ss_assemblage);
if (use.Get_ss_assemblage_ptr() != NULL)
{
size_t count_ss = use.Get_ss_assemblage_ptr()->Get_SSs().size();
last_model.count_ss_assemblage = count_ss;
last_model.ss_assemblage =
(const char **) PHRQ_malloc(count_ss * sizeof(char *));
if (last_model.ss_assemblage == NULL)
malloc_error();
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t j = 0; j < ss_ptrs.size(); j++)
{
last_model.ss_assemblage[i] = string_hsave(ss_ptrs[j]->Get_name().c_str());
}
}
else
{
last_model.count_ss_assemblage = 0;
last_model.ss_assemblage = NULL;
}
#ifdef SKIP
if (use.Get_ss_assemblage_ptr() != NULL)
{
last_model.count_ss_assemblage = use.Get_ss_assemblage_ptr()->count_s_s;
last_model.ss_assemblage =
@ -5883,6 +6231,7 @@ save_model(void)
last_model.count_ss_assemblage = 0;
last_model.ss_assemblage = NULL;
}
#endif
/*
* save list of phase pointers for pp_assemblage
*/
@ -6110,6 +6459,25 @@ check_same_model(void)
/*
* Check solid solutions
*/
if (use.Get_ss_assemblage_ptr() != NULL)
{
if (last_model.count_ss_assemblage != (int) use.Get_ss_assemblage_ptr()->Get_SSs().size())
return (FALSE);
std::vector<cxxSS *> ss_ptrs = use.Get_ss_assemblage_ptr()->Vectorize();
for (size_t i = 0; i < ss_ptrs.size(); i++)
{
if (last_model.ss_assemblage[i] != string_hsave(ss_ptrs[i]->Get_name().c_str()))
{
return (FALSE);
}
}
}
else
{
if (last_model.ss_assemblage != NULL)
return (FALSE);
}
#ifdef SKIP
if (use.Get_ss_assemblage_ptr() != NULL)
{
if (last_model.count_ss_assemblage !=
@ -6129,6 +6497,7 @@ check_same_model(void)
if (last_model.ss_assemblage != NULL)
return (FALSE);
}
#endif
/*
* Check pure_phases
*/

564
phreeqc/read.cpp
View File

@ -13,6 +13,9 @@
#include "GasPhase.h"
#include "Reaction.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "SS.h"
#include "SScomp.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
@ -9204,7 +9207,566 @@ read_user_graph(void)
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_solid_solutions(void)
/* ---------------------------------------------------------------------- */
{
/*
* Reads solid solution data
*
* Arguments:
* none
*
* Returns:
* KEYWORD if keyword encountered, input_error may be incremented if
* a keyword is encountered in an unexpected position
* EOF if eof encountered while reading mass balance concentrations
* ERROR if error occurred reading data
*
*/
//int i, j, n, l;
//int count_s_s, number_s_s, count_comps;
int n_user, n_user_end;
char *ptr;
char *description;
//char token[MAX_LENGTH];
std::string token;
int return_value, opt;
char *next_char;
const char *opt_list[] = {
"component", /* 0 */
"comp", /* 1 */
"parms", /* 2 */
"gugg_nondimensional", /* 3 */
"gugg_kj", /* 4 */
"activity_coefficients", /* 5 */
"distribution_coefficients", /* 6 */
"miscibility_gap", /* 7 */
"spinodal_gap", /* 8 */
"critical_point", /* 9 */
"alyotropic_point", /* 10 */
"temp", /* 11 */
"tempk", /* 12 */
"tempc", /* 13 */
"thompson", /* 14 */
"margules", /* 15 */
"comp1", /* 16 */
"comp2" /* 17 */
};
int count_opt_list = 18;
/*
* Read ss_assemblage number
*/
ptr = line;
read_number_description(ptr, &n_user, &n_user_end, &description);
cxxSSassemblage temp_ss_assemblage;
temp_ss_assemblage.Set_n_user(n_user);
temp_ss_assemblage.Set_n_user_end(n_user_end);
temp_ss_assemblage.Set_description(description);
free_check_null(description);
temp_ss_assemblage.Set_new_def(true);
std::vector<cxxSScomp> comps;
cxxSScomp * comp0_ptr = NULL;
cxxSScomp * comp1_ptr = NULL;
cxxSS * ss_ptr = NULL;
/*
* Set use data to first read
*/
if (!use.Get_ss_assemblage_in())
{
use.Set_ss_assemblage_in(true);
use.Set_n_ss_assemblage_user(n_user);
}
/*
* Read solid solutions
*/
return_value = UNKNOWN;
for (;;)
{
opt = get_option(opt_list, count_opt_list, &next_char);
switch (opt)
{
case OPTION_EOF: /* end of file */
return_value = EOF;
break;
case OPTION_KEYWORD: /* keyword */
return_value = KEYWORD;
break;
case OPTION_ERROR:
input_error++;
error_msg("Unknown input in SOLID_SOLUTIONS keyword.", CONTINUE);
error_msg(line_save, CONTINUE);
break;
/*
* New component
*/
case 0: /* component */
case 1: /* comp */
{
cxxSScomp comp;
/*
* Read phase name of component
*/
ptr = next_char;
copy_token(token, &ptr);
comp.Set_name(token);
/*
* Read moles of component
*/
if (copy_token(token, &ptr) == EMPTY)
{
comp.Set_moles(NAN);
}
else
{
int j = sscanf(token.c_str(), SCANFORMAT, &dummy);
comp.Set_moles(dummy);
(LDBLE) dummy;
if (j != 1)
{
error_msg("Expected moles of solid solution.", CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
}
}
comps.push_back(comp);
}
break;
case 2: /* parms */
case 3: /* gugg_nondimensional */
/*
* Read parameters
*/
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p()[0] = dummy;
}
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p()[1] = dummy;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_A0_A1);
break;
case 4: /* gugg_kj */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p()[0] = dummy;
}
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p()[1] = dummy;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_DIM_GUGG);
break;
case 5: /* activity coefficients */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 4; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 4)
{
error_string = sformatf(
"Expected 4 parameters to calculate a0 and a1 from two activity coefficients, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_GAMMAS);
break;
case 6: /* distribution coefficients */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 4; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 4)
{
error_string = sformatf(
"Expected 4 parameters to calculate a0 and a1 from two distribution coefficients, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_DIST_COEF);
break;
case 7: /* miscibility_gap */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected 2 miscibility gap fractions of component 2 to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_MISCIBILITY);
break;
case 8: /* spinodal_gap */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected 2 spinodal gap fractions of component 2 to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_SPINODAL);
break;
case 9: /* critical point */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected fraction of component 2 and critical temperature to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_CRITICAL);
break;
case 10: /* alyotropic point */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected fraction of component 2 and sigma pi at alyotropic point to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_ALYOTROPIC);
break;
case 12: /* tempk */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
{
ptr = next_char;
int j = 0;
if (copy_token(token, &ptr) != EMPTY)
{
j = sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Set_tk(dummy);
}
if (j != 1)
{
error_string = sformatf(
"Expected temperature (Kelvin) for parameters, assemblage %d, solid solution %s, using 298.15 K",
n_user,
ss_ptr->Get_name().c_str());
warning_msg(error_string);
ss_ptr->Set_tk(298.15);
}
}
break;
case 11: /* temp */
case 13: /* tempc */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
{
ptr = next_char;
int j = 0;
if (copy_token(token, &ptr) != EMPTY)
{
j = sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Set_tk(dummy + 298.15);
}
if (j != 1)
{
error_string = sformatf(
"Expected temperature (Celcius) for parameters, assemblage %d, solid solution %s, using 25 C",
n_user,
ss_ptr->Get_name().c_str());
warning_msg(error_string);
ss_ptr->Set_tk(298.15);
}
}
break;
case 14: /* Thompson and Waldbaum */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected Wg2 and Wg1 Thompson-Waldbaum parameters to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_WALDBAUM);
break;
case 15: /* Margules */
if (!ss_ptr)
{
error_msg("Solid solution name has not been defined", CONTINUE);
break;
}
ptr = next_char;
ss_ptr->Get_p().clear();
for (int i = 0; i < 2; i++)
{
if (copy_token(token, &ptr) != EMPTY)
{
sscanf(token.c_str(), SCANFORMAT, &dummy);
ss_ptr->Get_p().push_back(dummy);
}
}
if (ss_ptr->Get_p().size() != 2)
{
error_string = sformatf(
"Expected alpha2 and alpha3 Margules parameters to calculate a0 and a1, assemblage %d, solid solution %s",
n_user,
ss_ptr->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
ss_ptr->Set_input_case(cxxSS::SS_PARM_MARGULES);
break;
case 16: /* comp1 */
/*
* Read phase name of component
*/
delete comp0_ptr;
comp0_ptr = new cxxSScomp;
ptr = next_char;
copy_token(token, &ptr);
comp0_ptr->Set_name(token);
/*
* Read moles of component
*/
if (copy_token(token, &ptr) == EMPTY)
{
comp0_ptr->Set_moles(NAN);
}
else
{
int j = sscanf(token.c_str(), SCANFORMAT, &dummy);
comp0_ptr->Set_moles(dummy);
if (j != 1)
{
error_msg("Expected moles of solid solution.", CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
}
}
break;
case 17: /* comp2 */
delete comp1_ptr;
comp1_ptr = new cxxSScomp;
/*
* Read phase name of component
*/
ptr = next_char;
copy_token(token, &ptr);
comp1_ptr->Set_name(token);
/*
* Read moles of component
*/
if (copy_token(token, &ptr) == EMPTY)
{
comp1_ptr->Set_moles(NAN);
}
else
{
int j = sscanf(token.c_str(), SCANFORMAT, &dummy);
comp1_ptr->Set_moles(dummy);
if (j != 1)
{
error_msg("Expected moles of solid solution.", CONTINUE);
error_msg(line_save, CONTINUE);
input_error++;
}
}
break;
/*
* New solid solution
*/
case OPTION_DEFAULT:
if(ss_ptr)
{
comps.insert(comps.begin(), *comp1_ptr);
comps.insert(comps.begin(), *comp0_ptr);
ss_ptr->Set_ss_comps(comps);
temp_ss_assemblage.Get_SSs()[ss_ptr->Get_name()] = *ss_ptr;
delete ss_ptr;
ss_ptr = NULL;
comps.clear();
delete comp0_ptr, comp1_ptr;
comp0_ptr = comp1_ptr = NULL;
}
ss_ptr = new cxxSS;
/*
* Read solid solution name
*/
ptr = line;
copy_token(token, &ptr);
ss_ptr->Set_name(token);
ss_ptr->Set_total_moles(NAN);
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
}
// add last ss and clean up
comps.insert(comps.begin(), *comp1_ptr);
comps.insert(comps.begin(), *comp0_ptr);
ss_ptr->Set_ss_comps(comps);
temp_ss_assemblage.Get_SSs()[ss_ptr->Get_name()] = *ss_ptr;
delete ss_ptr;
ss_ptr = NULL;
comps.clear();
delete comp0_ptr, comp1_ptr;
comp0_ptr = comp1_ptr = NULL;
// check non ideal ss
std::vector<cxxSS *> ss_v;
for (size_t i = 0; i < ss_v.size(); i++)
{
if (ss_v[i]->Get_p()[0] != 0.0 ||
ss_v[i]->Get_p()[1] != 0.0)
{
if (ss_v[i]->Get_ss_comps().size() != 2)
{
error_string = sformatf(
"Solid solution, %s, is nonideal. Must define exactly two components (-comp1 and -comp2).",
ss_v[i]->Get_name().c_str());
error_msg(error_string, CONTINUE);
input_error++;
}
}
}
// Add to map
Rxn_ss_assemblage_map[n_user] = temp_ss_assemblage;
return (return_value);
}
#ifdef SKIP
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_solid_solutions(void)
@ -9763,7 +10325,7 @@ read_solid_solutions(void)
(size_t) count_s_s, (size_t) sizeof(struct s_s), s_s_compare);
return (return_value);
}
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
read_llnl_aqueous_model_parameters(void)

2
phreeqc/sit.cpp
View File

@ -845,7 +845,7 @@ jacobian_sit(void)
break;
case MU:
case PP:
case S_S_MOLES:
case SS_MOLES:
continue;
break;
}

51
phreeqc/structures.cpp
View File

@ -10,6 +10,7 @@
#include "Reaction.h"
#include "PPassemblage.h"
#include "Use.h"
#include "SSassemblage.h"
/* ---------------------------------------------------------------------- */
int Phreeqc::
@ -104,14 +105,15 @@ clean_up(void)
#endif
/* s_s assemblages */
Rxn_ss_assemblage_map.clear();
#ifdef SKIP
for (j = 0; j < count_ss_assemblage; j++)
{
ss_assemblage_free(&ss_assemblage[j]);
}
ss_assemblage =
(struct ss_assemblage *) free_check_null(ss_assemblage);
#endif
/* irreversible reactions */
Rxn_reaction_map.clear();
@ -375,7 +377,7 @@ clean_up(void)
count_surface = 0;
//count_gas_phase = 0;
count_kinetics = 0;
count_ss_assemblage = 0;
//count_ss_assemblage = 0;
count_elements = 0;
//count_irrev = 0;
@ -436,13 +438,14 @@ reinitialize(void)
count_pp_assemblage = 0;
#endif
/* s_s assemblages */
Rxn_ss_assemblage_map.clear();
#ifdef SKIP
for (j = 0; j < count_ss_assemblage; j++)
{
ss_assemblage_free(&ss_assemblage[j]);
}
count_ss_assemblage = 0;
#endif
/* gases */
Rxn_gas_phase_map.clear();
@ -3233,7 +3236,7 @@ s_store(const char *name, LDBLE l_z, int replace_if_found)
return (s_ptr);
}
#ifdef SKIP
/* **********************************************************************
*
* Routines related to structure "ss_assemblage"
@ -3656,7 +3659,7 @@ s_s_compare(const void *ptr1, const void *ptr2)
return (strcmp_nocase(s_s_ptr1->name, s_s_ptr2->name));
}
#endif
/* **********************************************************************
*
* Routines related to structure "save_values"
@ -6451,7 +6454,7 @@ cxxSolutionIsotopeList2isotope(const cxxSolutionIsotopeList * il)
}
return (iso);
}
#ifdef SKIP
#include "../SSassemblage.h"
#include "../SS.h"
struct ss_assemblage * Phreeqc::
@ -6554,7 +6557,7 @@ cxxSSassemblageSS2s_s(const std::map < std::string, cxxSS > * sscomp)
}
return (s_s_ptr);
}
#endif
#include "../Surface.h"
struct surface * Phreeqc::
cxxSurface2surface(const cxxSurface * surf)
@ -6958,6 +6961,15 @@ Use2cxxStorageBin(cxxStorageBin & sb)
sb.Set_GasPhase(use.Get_n_gas_phase_user(), entity_ptr);
}
}
if (use.Get_ss_assemblage_in())
{
cxxSSassemblage *entity_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, use.Get_n_ss_assemblage_user());
if (entity_ptr != NULL)
{
sb.Set_SSassemblage(use.Get_n_ss_assemblage_user(), entity_ptr);
}
}
#ifdef SKIP
if (use.Get_ss_assemblage_in())
{
struct ss_assemblage *struct_entity = ss_assemblage_bsearch(use.Get_n_ss_assemblage_user(), &n);
@ -6967,6 +6979,7 @@ Use2cxxStorageBin(cxxStorageBin & sb)
sb.Set_SSassemblage(use.Get_n_ss_assemblage_user(), &entity);
}
}
#endif
if (use.Get_kinetics_in())
{
struct kinetics *struct_entity = kinetics_bsearch(use.Get_n_kinetics_user(), &n);
@ -7059,12 +7072,20 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb)
}
#endif
// SSassemblages
{
std::map<int, cxxSSassemblage>::iterator it;
for (it = Rxn_ss_assemblage_map.begin(); it != Rxn_ss_assemblage_map.end(); it++)
{
sb.Set_SSassemblage(it->second.Get_n_user(), &(it->second));
}
}
#ifdef SKIP
for (i = 0; i < count_ss_assemblage; i++)
{
cxxSSassemblage entity(&ss_assemblage[i], sb.Get_io());
sb.Set_SSassemblage(ss_assemblage[i].n_user, &entity );
}
#endif
// Surfaces
for (i = 0; i < count_surface; i++)
{
@ -7172,6 +7193,14 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb, int n)
}
#endif
// SSassemblages
{
cxxSSassemblage *entity_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, n);
if (entity_ptr != NULL)
{
sb.Set_SSassemblage(n, entity_ptr);
}
}
#ifdef SKIP
{
if (ss_assemblage_bsearch(n, &pos) != NULL)
{
@ -7180,7 +7209,7 @@ phreeqc2cxxStorageBin(cxxStorageBin & sb, int n)
sb.Set_SSassemblage(n, &ent);
}
}
#endif
// Surfaces
{
if (surface_bsearch(n, &pos) != NULL)