Working on edl for PHREEQCRM

Adding edl_species for PHREEQC.

git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/trunk@9915 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2015-06-30 17:06:53 +00:00
parent f15d9155fc
commit 49ba8bca71
5 changed files with 245 additions and 3 deletions

View File

@ -1445,6 +1445,10 @@ listtokens(FILE * f, tokenrec * l_buf)
output_msg("SURF");
break;
case tokedl_species:
output_msg("EDL_SPECIES");
break;
case tokstep_no:
output_msg("STEP_NO");
break;
@ -2677,6 +2681,107 @@ factor(struct LOC_exec * LINK)
}
break;
case tokedl_species:
{
double area, thickness;
require(toklp, LINK);
const char *surf_name = stringfactor(STR1, LINK);
require(tokcomma, LINK);
// variable for number of species
LINK->t = LINK->t->next;
count_varrec = LINK->t->UU.vp;
if (LINK->t->kind != tokvar || !count_varrec || count_varrec->stringvar != 0)
{
snerr(": can`t find variable");
exit(4);
}
require(tokcomma, LINK);
// variable for species names
LINK->t = LINK->t->next;
require(tokcomma, LINK);
names_varrec = LINK->t->UU.vp;
if (LINK->t->kind != tokvar || !names_varrec || names_varrec->stringvar != 1)
{
snerr(": can`t find name of species");
exit(4);
}
// variable for species concentrations
LINK->t = LINK->t->next;
require(tokcomma, LINK);
moles_varrec = LINK->t->UU.vp;
if (LINK->t->kind != tokvar || moles_varrec->stringvar != 0)
snerr(": can`t find concentrations of species");
LINK->t = LINK->t->next;
// variable for area
LINK->t = LINK->t->next;
require(tokcomma, LINK);
varrec *area_varrec = LINK->t->UU.vp;
if (LINK->t->kind != tokvar || area_varrec->stringvar != 0)
snerr(": can`t find area varaiable");
LINK->t = LINK->t->next;
// varaiable for thickness
LINK->t = LINK->t->next;
require(tokcomma, LINK);
varrec *thickness_varrec = LINK->t->UU.vp;
if (LINK->t->kind != tokvar || thickness_varrec->stringvar != 0)
snerr(": can`t find thickness variable");
LINK->t = LINK->t->next;
require(tokrp, LINK);
free_dim_stringvar(names_varrec);
PhreeqcPtr->free_check_null(moles_varrec->UU.U0.arr);
moles_varrec->UU.U0.arr = NULL;
// Call subroutine
if (parse_all)
{
PhreeqcPtr->sys_tot = 0;
PhreeqcPtr->count_sys = 1000;
int count_sys = PhreeqcPtr->count_sys;
names_arg = (char **) PhreeqcPtr->PHRQ_calloc((size_t) (count_sys + 1), sizeof(char *));
if (names_arg == NULL)
{
PhreeqcPtr->malloc_error();
exit(4);
}
moles_arg = (LDBLE *) PhreeqcPtr->PHRQ_calloc((size_t) (count_sys + 1), sizeof(LDBLE));
if (moles_arg == NULL)
{
PhreeqcPtr->malloc_error();
exit(4);
}
names_arg[0] = NULL;
moles_arg[0] = 0;
count_species = (LDBLE) count_sys;
n.UU.val = 0;
}
else
{
//n.UU.val = PhreeqcPtr->system_total(elt_name, &count_species, &(names_arg),
// &(types_arg), &(moles_arg));
n.UU.val = PhreeqcPtr->edl_species(surf_name, &count_species, &(names_arg), &(moles_arg), &area, &thickness);
}
/*
* fill in varrec structures
*/
*count_varrec->UU.U0.val = count_species;
names_varrec->UU.U1.sarr = names_arg;
moles_varrec->UU.U0.arr = moles_arg;
*area_varrec->UU.U0.val = area;
*thickness_varrec->UU.U0.val = thickness;
for (i = 0; i < maxdims; i++)
{
names_varrec->dims[i] = 0;
moles_varrec->dims[i] = 0;
}
names_varrec->dims[0] = (long) (*count_varrec->UU.U0.val) + 1;
moles_varrec->dims[0] = (long) (*count_varrec->UU.U0.val) + 1;
names_varrec->numdims = 1;
moles_varrec->numdims = 1;
}
break;
case toklist_s_s:
{
/* list_s_s("calcite", count, name$, moles) */
@ -6970,7 +7075,8 @@ const std::map<const std::string, PBasic::BASIC_TOKEN>::value_type temp_tokens[]
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("equiv_frac", PBasic::tokeq_frac),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("callback", PBasic::tokcallback),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("diff_c", PBasic::tokdiff_c),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sa_declercq", PBasic::toksa_declercq)
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("sa_declercq", PBasic::toksa_declercq),
std::map<const std::string, PBasic::BASIC_TOKEN>::value_type("edl_species", PBasic::tokedl_species)
};
std::map<const std::string, PBasic::BASIC_TOKEN> PBasic::command_tokens(temp_tokens, temp_tokens + sizeof temp_tokens / sizeof temp_tokens[0]);

View File

@ -323,7 +323,8 @@ public:
tokequiv_frac,
tokcallback,
tokdiff_c,
toksa_declercq
toksa_declercq,
tokedl_species
};
#if !defined(PHREEQCI_GUI)

View File

@ -111,6 +111,8 @@ public:
LDBLE calc_logk_s(const char *name);
LDBLE calc_surface_charge(const char *surface_name);
LDBLE diff_layer_total(const char *total_name, const char *surface_name);
LDBLE edl_species(const char *surf_name, LDBLE * count, char ***names, LDBLE ** moles, LDBLE * area, LDBLE * thickness);
int get_edl_species(cxxSurfaceCharge & charge_ref);
LDBLE equi_phase(const char *phase_name);
LDBLE equi_phase_delta(const char *phase_name);
LDBLE equivalent_fraction(const char *name, LDBLE *eq, std::string &elt_name);

View File

@ -2581,7 +2581,133 @@ total_mole(const char *total_name)
}
return (t);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
get_edl_species(cxxSurfaceCharge & charge_ref)
/* ---------------------------------------------------------------------- */
{
double mass_water_surface = charge_ref.Get_mass_water();
space((void **) ((void *) &sys), count_s_x, &max_sys, sizeof(struct system_species));
count_sys = 0;
for (int j = 0; j < count_s_x; j++)
{
if (s_x[j]->type == H2O)
{
sys[count_sys].name = string_duplicate(s_x[j]->name);
sys[count_sys].moles = mass_water_surface / gfw_water;
sys_tot += sys[count_sys].moles;
}
else if (s_x[j]->type < H2O)
{
double molality = under(s_x[j]->lm);
double moles_excess = mass_water_aq_x * molality * charge_ref.Get_g_map()[s_x[j]->z].Get_g();
double moles_surface = mass_water_surface * molality + moles_excess;
sys[count_sys].name = string_duplicate(s_x[j]->name);
sys[count_sys].moles = moles_surface;
sys_tot += sys[count_sys].moles;
count_sys++;
#ifdef SKIP
double g = charge_ref.Get_g_map()[s_x[j]->z].Get_g();
double moles_excess = mass_water_aq_x * molality * (g * s_x[j]->erm_ddl +
mass_water_surface /
mass_water_aq_x * (s_x[j]->erm_ddl - 1));
double c = (mass_water_surface * molality + moles_excess) / mass_water_surface;
charge_ref.Get_dl_species_map()[s_x[j]->number] = c;
#endif
}
else
{
continue;
}
}
#ifdef SKIP
/*
* Provides total moles in system and lists of species/phases in sort order
*/
int i;
/*
* find total moles in aq, surface, and exchange
*/
for (i = 0; i < count_s_x; i++)
{
//if (s_x[i]->type != AQ)
if (s_x[i]->type > HPLUS)
continue;
sys[count_sys].name = string_duplicate(s_x[i]->name);
sys[count_sys].moles = s_x[i]->moles;
sys_tot += sys[count_sys].moles;
sys[count_sys].type = string_duplicate("aq");
count_sys++;
space((void **) ((void *) &sys), count_sys, &max_sys,
sizeof(struct system_species));
}
#endif
return (OK);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
edl_species(const char *surf_name, LDBLE * count, char ***names, LDBLE ** moles, LDBLE * area, LDBLE * thickness)
/* ---------------------------------------------------------------------- */
{
/*
* Provides total moles in system and lists of species/phases in sort order
*/
int i;
sys_tot = 0;
count_sys = 0;
max_sys = 100;
space((void **) ((void *) &sys), INIT, &max_sys,
sizeof(struct system_species));
if (!(dl_type_x == cxxSurface::NO_DL))
{
cxxSurface *surface_ptr = use.Get_surface_ptr();
for (size_t i = 0; i < surface_ptr->Get_surface_charges().size(); i++)
{
cxxSurfaceCharge & charge_ref = surface_ptr->Get_surface_charges()[i];
if (strcmp(charge_ref.Get_name().c_str(), surf_name) == 0)
{
get_edl_species(charge_ref);
*area = charge_ref.Get_specific_area() * charge_ref.Get_grams();
*thickness = surface_ptr->Get_thickness();
break;
}
}
}
/*
* Sort system species
*/
if (count_sys > 1)
{
qsort(sys, (size_t) count_sys,
(size_t) sizeof(struct system_species), system_species_compare);
}
/*
* malloc space
*/
*names = (char **) PHRQ_malloc((size_t) (count_sys + 1) * sizeof(char *));
if (names == NULL)
malloc_error();
*moles = (LDBLE *) PHRQ_malloc((size_t) (count_sys + 1) * sizeof(LDBLE));
if (moles == NULL)
malloc_error();
(*names)[0] = NULL;
(*moles)[0] = 0;
for (i = 0; i < count_sys; i++)
{
(*names)[i + 1] = sys[i].name;
(*moles)[i + 1] = sys[i].moles;
}
*count = (LDBLE) count_sys;
PHRQ_free(sys);
return (sys_tot);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
system_total(const char *total_name, LDBLE * count, char ***names,
char ***types, LDBLE ** moles)

View File

@ -1749,16 +1749,23 @@ xsurface_save(int n_user)
double mass_water_surface = charge_ref.Get_mass_water();
for (int j = 0; j < count_s_x; j++)
{
if (s_x[j]->type > HPLUS)
if (s_x[j]->type > H2O)
continue;
double molality = under(s_x[j]->lm);
double moles_excess = mass_water_aq_x * molality * charge_ref.Get_g_map()[s_x[j]->z].Get_g();
double moles_surface = mass_water_surface * molality + moles_excess;
charge_ref.Get_dl_species_map()[s_x[j]->number] = moles_surface/mass_water_surface;
#ifdef SKIP
double g = charge_ref.Get_g_map()[s_x[j]->z].Get_g();
double moles_excess = mass_water_aq_x * molality * (g * s_x[j]->erm_ddl +
mass_water_surface /
mass_water_aq_x * (s_x[j]->erm_ddl - 1));
double c = (mass_water_surface * molality + moles_excess) / mass_water_surface;
charge_ref.Get_dl_species_map()[s_x[j]->number] = c;
#endif
}
//charge_ref.Get_dl_species_map()[s_h2o->number] = 0.0;
charge_ref.Get_dl_species_map()[s_h2o->number] = 1.0/gfw_water;
}
}