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 '6a49d41253c7705c0273262d0290389b18643d32'

Browse Source
This commit is contained in:
Scott R Charlton 2019-10-29 12:29:32 -06:00
commit 63108ac76d
25 changed files with 1062 additions and 390 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
phreeqcpp/ExchComp.cxx
View File

@ -360,7 +360,8 @@ cxxExchComp::add(const cxxExchComp & addee, LDBLE extensive)
this->la = f1 * this->la + f2 * addee.la;
this->charge_balance += addee.charge_balance * extensive;
if (this->phase_name != addee.phase_name)
//if (this->phase_name != addee.phase_name)
if (Phreeqc::strcmp_nocase(this->phase_name.c_str(), addee.phase_name.c_str()) != 0)
{
std::ostringstream oss;
oss <<
@ -374,7 +375,8 @@ cxxExchComp::add(const cxxExchComp & addee, LDBLE extensive)
this->phase_proportion =
this->phase_proportion * f1 + addee.phase_proportion * f2;
}
if (this->rate_name != addee.rate_name)
//if (this->rate_name != addee.rate_name)
if (Phreeqc::strcmp_nocase(this->rate_name.c_str(), addee.rate_name.c_str()) != 0)
{
std::ostringstream oss;
oss <<

2
phreeqcpp/Exchange.cxx
View File

@ -467,7 +467,7 @@ cxxExchange::Deserialize(Dictionary & dictionary, std::vector < int >&ints, std:
this->exchange_comps.clear();
for (int n = 0; n < count; n++)
{
cxxExchComp ec;
cxxExchComp ec(this->io);
ec.Deserialize(dictionary, ints, doubles, ii, dd);
this->exchange_comps.push_back(ec);
}

2
phreeqcpp/PPassemblage.cxx
View File

@ -373,7 +373,7 @@ cxxPPassemblage::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->pp_assemblage_comps.clear();
for (int n = 0; n < count; n++)
{
cxxPPassemblageComp ppc;
cxxPPassemblageComp ppc(this->io);
ppc.Deserialize(dictionary, ints, doubles, ii, dd);
std::string str(ppc.Get_name());
this->pp_assemblage_comps[str] = ppc;

2
phreeqcpp/Phreeqc.cpp
View File

@ -713,6 +713,7 @@ void Phreeqc::init(void)
interlayer_Dflag = FALSE;
implicit = FALSE;
max_mixf = 1.0;
min_dif_LM = -30.0;
default_Dw = 0;
correct_Dw = 0;
multi_Dpor = 0;
@ -1691,6 +1692,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
interlayer_Dflag = pSrc->interlayer_Dflag;
implicit = pSrc->implicit;
max_mixf = pSrc->max_mixf;
min_dif_LM = pSrc->min_dif_LM;
default_Dw = pSrc->default_Dw;
correct_Dw = pSrc->correct_Dw;
multi_Dpor = pSrc->multi_Dpor;

9
phreeqcpp/Phreeqc.h
View File

@ -1066,9 +1066,11 @@ public:
LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant);
void diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant);
int fill_spec(int cell_no);
void define_ct_structures(void);
void diffuse_implicit(LDBLE DDt, int stagnant);
int fill_spec(int cell_no, int ref_cell);
LDBLE moles_from_redox_states(cxxSolution *sptr, const char *name);
LDBLE moles_from_donnan_layer(cxxSurface *sptr, const char *name, LDBLE moles_needed);
LDBLE add_MCD_moles(LDBLE moles, LDBLE min_mol, int i, cxxSolution *sptr, const char *name);
int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i, int stagnant);
static int sort_species_name(const void *ptr1, const void *ptr2);
int disp_surf(LDBLE stagkin_time);
@ -1429,6 +1431,7 @@ protected:
int interlayer_Dflag; /* multicomponent diffusion and diffusion through interlayer porosity */
int implicit; /* implicit calculation of diffusion */
LDBLE max_mixf; /* the maximum value of the implicit mixfactor = De * Dt / (Dx^2) */
LDBLE min_dif_LM; /* the minimal log10(molality) for including a species in multicomponent diffusion */
LDBLE default_Dw; /* default species diffusion coefficient in water at 25oC, m2/s */
int correct_Dw; /* if true, Dw is adapted in calc_SC */
LDBLE multi_Dpor; /* uniform porosity of free porewater in solid medium */

2
phreeqcpp/SS.cxx
View File

@ -580,7 +580,7 @@ cxxSS::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->ss_comps.clear();
for (int i = 0; i < count; i++)
{
cxxSScomp ssc;
cxxSScomp ssc(this->io);
ssc.Deserialize(dictionary, ints, doubles, ii, dd);
this->ss_comps.push_back(ssc);
}

2
phreeqcpp/SSassemblage.cxx
View File

@ -319,7 +319,7 @@ cxxSSassemblage::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->SSs.clear();
for (int n = 0; n < count; n++)
{
cxxSS ssc;
cxxSS ssc(this->io);
ssc.Deserialize(dictionary, ints, doubles, ii, dd);
std::string str(ssc.Get_name());
this->SSs[str] = ssc;

6
phreeqcpp/Serializer.cxx
View File

@ -133,7 +133,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break;
case PT_GASPHASE:
{
cxxGasPhase entity;
cxxGasPhase entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_gas_phase_map()[n_user] = entity;
@ -141,7 +141,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break;
case PT_KINETICS:
{
cxxKinetics entity;
cxxKinetics entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_kinetics_map()[n_user] = entity;
@ -158,7 +158,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break;
case PT_SSASSEMBLAGE:
{
cxxSSassemblage entity;
cxxSSassemblage entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_ss_assemblage_map()[n_user] = entity;

4
phreeqcpp/StorageBin.cxx
View File

@ -1061,7 +1061,7 @@ cxxStorageBin::read_raw(CParser & parser)
break;
case Keywords::KEY_SOLID_SOLUTIONS_RAW:
{
cxxSSassemblage entity;
cxxSSassemblage entity(this->Get_io());
entity.read_raw(parser);
SSassemblages[entity.Get_n_user()] = entity;
}
@ -1195,7 +1195,7 @@ cxxStorageBin::read_raw_keyword(CParser & parser)
case Keywords::KEY_SOLID_SOLUTIONS_RAW:
{
cxxSSassemblage entity;
cxxSSassemblage entity(this->Get_io());
entity.read_raw(parser);
SSassemblages[entity.Get_n_user()] = entity;
entity_number = entity.Get_n_user();

4
phreeqcpp/Surface.cxx
View File

@ -812,7 +812,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->surface_comps.clear();
for (int n = 0; n < count; n++)
{
cxxSurfaceComp sc;
cxxSurfaceComp sc(this->io);
sc.Deserialize(dictionary, ints, doubles, ii, dd);
this->surface_comps.push_back(sc);
}
@ -822,7 +822,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->surface_charges.clear();
for (int n = 0; n < count; n++)
{
cxxSurfaceCharge sc;
cxxSurfaceCharge sc(this->io);
sc.Deserialize(dictionary, ints, doubles, ii, dd);
this->surface_charges.push_back(sc);
}

10
phreeqcpp/SurfaceComp.cxx
View File

@ -395,7 +395,8 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
this->totals.add_extensive(addee.totals, extensive);
this->la = f1 * this->la + f2 * addee.la;
this->charge_balance += addee.charge_balance * extensive;
if (this->phase_name != addee.phase_name)
//if (this->phase_name != addee.phase_name)
if (Phreeqc::strcmp_nocase(this->phase_name.c_str(),addee.phase_name.c_str()) != 0)
{
std::ostringstream oss;
oss <<
@ -410,11 +411,12 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
this->phase_proportion * f1 + addee.phase_proportion * f2;
}
if (this->rate_name != addee.rate_name)
//if (this->rate_name != addee.rate_name)
if (Phreeqc::strcmp_nocase(this->rate_name.c_str(), addee.rate_name.c_str()) != 0)
{
std::ostringstream oss;
oss <<
"Cannot mix two exchange components with same formula and different related kinetics, "
"Cannot mix two surface components with same formula and different related kinetics, "
<< this->formula;
error_msg(oss.str().c_str(), CONTINUE);
return;
@ -430,7 +432,7 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
{
std::ostringstream oss;
oss <<
"Cannot mix exchange components related to phase with exchange components related to kinetics, "
"Cannot mix surface components related to phase with exchange components related to kinetics, "
<< this->formula;
error_msg(oss.str().c_str(), CONTINUE);
return;

10
phreeqcpp/class_main.cpp
View File

@ -277,8 +277,7 @@ write_banner(void)
/* version */
#ifdef NPP
//len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1");
len = sprintf(buffer, "* PHREEQC-%s *", "3.5.2");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -302,7 +301,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "May 29, 2019");
len = sprintf(buffer, "%s", "August 1, 2019");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -491,13 +490,10 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
user_database = string_duplicate(token);
screen_msg(sformatf("Database file: %s\n\n", token));
strcpy(db_file, token);
#ifdef NPP
//output_msg(sformatf("Using PHREEQC: version 3.5.1 Ampère, compiled May 29, 2019\n"));
#endif
output_msg(sformatf(" Input file: %s\n", in_file));
output_msg(sformatf(" Output file: %s\n", out_file));
#ifdef NPP
output_msg(sformatf("Using PHREEQC: version 3.5.1, compiled May 29, 2019\n"));
output_msg(sformatf("Using PHREEQC: version 3.5.2, compiled August 1, 2019\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token));
#ifdef NPP

1
phreeqcpp/common/PHRQ_io.cpp
View File

@ -792,6 +792,7 @@ get_line(void)
continue;
#else
output_msg(errstr.str().c_str());
output_flush();
error_msg(errstr.str().c_str(), OT_STOP);
#endif
}

81
phreeqcpp/cvode.cpp
View File

@ -72,6 +72,7 @@
#include "sundialsmath.h"
#include "Phreeqc.h"
#if !defined(WIN32_MEMORY_DEBUG)
#define malloc MACHENV_MALLOC PHRQ_malloc
#endif
@ -1974,7 +1975,6 @@ CVStep(CVodeMem cv_mem)
}
#endif
CVSet(cv_mem);
nflag = CVnls(cv_mem, nflag);
if (CVMEM cvode_error == TRUE || predict_fail)
{
@ -2785,7 +2785,6 @@ CVnlsNewton(CVodeMem cv_mem, int nflag)
loop
{
f(N, tn, zn[0], ftemp, f_data);
nfe++;
@ -3408,7 +3407,7 @@ CVHandleFailure(CVodeMem cv_mem, int kflag)
size is reset accordingly.
*****************************************************************/
#ifdef ORIGINAL
void
CVBDFStab(CVodeMem cv_mem)
{
@ -3469,6 +3468,82 @@ CVBDFStab(CVodeMem cv_mem)
nscon = 0;
}
}
#endif
void
CVBDFStab(CVodeMem cv_mem)
{
// appt try...
if (q >= 3 && qprime >= q)
{
if (tq[5] < saved_tq5)
qprime = 1;
//else
//nscon = 0;
}
//int i, k, ldflag, factorial;
//realtype sq, sqm1, sqm2;
///* If order is 3 or greater, then save scaled derivative data,
// push old data down in i, then add current values to top. */
//if (q >= 3)
//{
// for (k = 1; k <= 3; k++)
// {
// for (i = 5; i >= 2; i--)
// ssdat[i][k] = ssdat[i - 1][k];
// }
// factorial = 1;
// for (i = 1; i <= q - 1; i++)
// factorial *= i;
// sq = factorial * q * (q + 1) * acnrm / MAX(tq[5], TINY);
// sqm1 = factorial * q * N_VWrmsNorm(zn[q], ewt);
// sqm2 = factorial * N_VWrmsNorm(zn[q - 1], ewt);
// ssdat[1][1] = sqm2 * sqm2;
// ssdat[1][2] = sqm1 * sqm1;
// ssdat[1][3] = sq * sq;
//}
//if (qprime >= q)
//{
// /* If order is 3 or greater, and enough ssdat has been saved,
// nscon >= q+5, then call stability limit detection routine. */
// if ((q >= 3) && (nscon >= q + 5))
// {
// ldflag = CVsldet(cv_mem);
// //cv_mem->cv_machenv->phreeqc_ptr->set_forward_output_to_log(1); // appt
// qprime = 1; // appt try
// //CVMEM warning_msg(CVMEM sformatf(
// // "CVBDFStab: ldflag = %d, order(q) = %d, qprime = %d, nst = %d, h = %8.2e, time = %8.2e\n",
// // ldflag, q, qprime, nst, h, CVMEM cvode_last_good_time));
// if (ldflag > 3)
// {
// /* A stability limit violation is indicated by
// a return flag of 4, 5, or 6.
// Reduce new order. */
// qprime = q - 1;
// eta = etaqm1;
// eta = MIN(eta, etamax);
// eta = eta / MAX(ONE, ABS(h) * hmax_inv * eta);
// hprime = h * eta;
// iopt[NOR] = iopt[NOR] + 1;
// //CVMEM warning_msg(CVMEM sformatf(
// // " Order reduced to %d by CVBDFStab at nst = %d,\n h = %e hnew = %e\n",
// // qprime,nst,h,h*eta));
// }
// }
//}
//else
//{
// /* Otherwise, let order increase happen, and
// reset stability limit counter, nscon. */
// nscon = 0;
//}
}
/********************* CVsldet ************************************
This routine detects stability limitation using stored scaled

2
phreeqcpp/cxxKinetics.cxx
View File

@ -605,7 +605,7 @@ cxxKinetics::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->kinetics_comps.clear();
for (int i = 0; i < n; i++)
{
cxxKineticsComp kc;
cxxKineticsComp kc(this->io);
kc.Deserialize(dictionary, ints, doubles, ii, dd);
this->kinetics_comps.push_back(kc);
}

5
phreeqcpp/global_structures.h
View File

@ -547,6 +547,7 @@ struct cell_data
LDBLE potV; /* potential (V) */
int punch;
int print;
int same_model;
};
/*----------------------------------------------------------------------
@ -1052,12 +1053,12 @@ struct sol_D
struct J_ij
{
const char *name;
LDBLE tot1, tot2, charge; /* species change in cells i and j */
LDBLE tot1, tot2, tot_stag, charge; /* species change in cells i and j */
};
struct M_S
{
const char *name;
LDBLE tot1, tot2, charge; /* master species transport in cells i and j */
LDBLE tot1, tot2, tot_stag, charge; /* master species transport in cells i and j */
};
// Pitzer definitions
typedef enum

12
phreeqcpp/integrate.cpp
View File

@ -1069,9 +1069,17 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
l_iter++;
if (l_iter > 50)
{
pr.all = TRUE;
pr.exchange = TRUE;
pr.headings = TRUE;
pr.pp_assemblage = TRUE;
pr.species = TRUE;
pr.surface = TRUE;
pr.totals = TRUE;
print_all();
error_string = sformatf(
"\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n",
(double) surf_chrg_eq, (double) charge_ptr->Get_mass_water());
"\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n",
(double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
error_msg(error_string, STOP);
}
}

10
phreeqcpp/kinetics.cpp
View File

@ -500,8 +500,8 @@ rk_kinetics(int i, LDBLE kin_time, int use_mix, int nsaver,
if (step_bad > kinetics_ptr->Get_bad_step_max())
{
error_string = sformatf(
"Bad RK steps > %d. Please decrease (time)step or increase -bad_step_max.",
kinetics_ptr->Get_bad_step_max());
"Bad RK steps > %d in cell %d. Please decrease (time)step or increase -bad_step_max.",
kinetics_ptr->Get_bad_step_max(), cell_no);
error_msg(error_string, STOP);
}
@ -1511,11 +1511,11 @@ restart:
}
if (j == 14)
{
cxxStorageBin error_bin;
cxxStorageBin error_bin(this->Get_phrq_io());
Use2cxxStorageBin(error_bin);
std::ostringstream error_input;
error_bin.dump_raw(error_input, 0);
cxxStorageBin reread;
cxxStorageBin reread(this->Get_phrq_io());
std::istringstream is(error_input.str());
CParser cp(is);
cp.set_echo_stream(CParser::EO_NONE);
@ -1574,7 +1574,7 @@ restart:
* write to error.inp what failed to converge.
*/
std::ofstream error_input("error.inp");
cxxStorageBin error_bin;
cxxStorageBin error_bin(this->Get_phrq_io());
Use2cxxStorageBin(error_bin);
error_bin.dump_raw(error_input, 0);
error_input.close();

6
phreeqcpp/mainsubs.cpp
View File

@ -1396,7 +1396,7 @@ xss_assemblage_save(int n_user)
* Save ss_assemblage composition into structure ss_assemblage with user
* number n_user.
*/
cxxSSassemblage temp_ss_assemblage;
cxxSSassemblage temp_ss_assemblage(this->phrq_io);
if (use.Get_ss_assemblage_ptr() == NULL)
return (OK);
@ -1719,6 +1719,8 @@ xsurface_save(int n_user)
if (x[i]->type == SURFACE)
{
cxxSurfaceComp *comp_ptr = temp_surface.Find_comp(x[i]->surface_comp);
if (comp_ptr == NULL)
continue; // appt in transport with different mobile and stagnant surfaces
assert(comp_ptr);
comp_ptr->Set_la(x[i]->master[0]->s->la);
comp_ptr->Set_moles(0.);
@ -1748,6 +1750,8 @@ xsurface_save(int n_user)
else if (x[i]->type == SURFACE_CB && (use.Get_surface_ptr()->Get_type() == cxxSurface::DDL || use.Get_surface_ptr()->Get_type() == cxxSurface::CCM))
{
cxxSurfaceCharge *charge_ptr = temp_surface.Find_charge(x[i]->surface_charge);
if (charge_ptr == NULL)
continue; // appt in transport with different mobile and stagnant surfaces
assert(charge_ptr);
charge_ptr->Set_charge_balance(x[i]->f);
charge_ptr->Set_la_psi(x[i]->master[0]->s->la);

4
phreeqcpp/phqalloc.cpp
View File

@ -141,8 +141,8 @@ PHRQ_calloc(size_t num, size_t size
assert((s_pTail == NULL) || (s_pTail->pNext == NULL));
p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num);
//p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
//p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num);
p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
if (p == NULL)
return NULL;

2
phreeqcpp/prep.cpp
View File

@ -6034,6 +6034,8 @@ check_same_model(void)
last_model.force_prep = FALSE;
return (FALSE);
}
if (state == TRANSPORT && cell_data[cell_no].same_model)
return TRUE;
/*
* Check master species
*/

5
phreeqcpp/print.cpp
View File

@ -1353,8 +1353,9 @@ print_pp_assemblage(void)
{
if (x[j]->type != PP)
continue;
//cxxPPassemblageComp * comp_ptr = pp_assemblage_ptr->Find(x[j]->pp_assemblage_comp_name);
cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr;
cxxPPassemblage * pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, use.Get_n_pp_assemblage_user());
cxxPPassemblageComp * comp_ptr = pp_assemblage_ptr->Find(x[j]->pp_assemblage_comp_name);
//cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr; // appt, is sometimes lost??
/*
* Print saturation index
*/

13
phreeqcpp/read.cpp
View File

@ -1007,7 +1007,7 @@ read_exchange(void)
input_error++;
break;
}
cxxExchComp temp_comp;
cxxExchComp temp_comp(this->phrq_io);
temp_exchange.Get_exchange_comps().push_back(temp_comp);
comp_ptr = &(temp_exchange.Get_exchange_comps().back());
comp_ptr->Set_formula(token.c_str());
@ -1268,7 +1268,7 @@ read_gas_phase(void)
char *ptr;
char *description;
char token[MAX_LENGTH];
cxxGasPhase temp_gas_phase;
cxxGasPhase temp_gas_phase(this->phrq_io);
int return_value, opt;
char *next_char;
const char *opt_list[] = {
@ -2057,7 +2057,7 @@ read_kinetics(void)
*/
ptr = line;
read_number_description(ptr, &n_user, &n_user_end, &description);
cxxKinetics temp_kinetics;
cxxKinetics temp_kinetics(this->phrq_io);
temp_kinetics.Set_n_user(n_user);
temp_kinetics.Set_n_user_end(n_user_end);
temp_kinetics.Set_description(description);
@ -7538,7 +7538,8 @@ read_surface(void)
break;
}
cxxSurfaceComp temp_comp;
cxxSurfaceComp temp_comp(this->phrq_io);
temp_surface.Get_surface_comps().push_back(temp_comp);
comp_ptr = &(temp_surface.Get_surface_comps().back());
comp_ptr->Set_formula(token.c_str());
@ -7648,7 +7649,7 @@ read_surface(void)
formula = (char*)free_check_null(formula);
if (charge_ptr == NULL)
{
cxxSurfaceCharge temp_charge;
cxxSurfaceCharge temp_charge(this->phrq_io);
temp_charge.Set_name(name);
if (comp_ptr->Get_phase_name().size() == 0
&& comp_ptr->Get_rate_name().size() == 0)
@ -10125,7 +10126,7 @@ read_solid_solutions(void)
case 0: /* component */
case 1: /* comp */
{
cxxSScomp comp;
cxxSScomp comp(this->phrq_io);
/*
* Read phase name of component
*/

63
phreeqcpp/readtr.cpp
View File

@ -37,13 +37,13 @@ read_transport(void)
*/
char *ptr;
int i, j, l;
int count_length, count_disp, count_punch, count_print, count_por;
int count_length, count_disp, count_punch, count_print, count_por, count_same_model;
int count_length_alloc, count_disp_alloc, count_por_alloc;
char token[MAX_LENGTH];
char *description;
int n_user, n_user_end;
LDBLE *length, *disp, *pors;
int *punch_temp, *print_temp;
int *punch_temp, *print_temp, *same_model_temp;
int return_value, opt, opt_save;
char *next_char, *next_char_save;
char file_name[MAX_LENGTH];
@ -95,9 +95,10 @@ read_transport(void)
"porosity", /* 43 */
"fix_current", /* 44 */
"current", /* 45 */
"implicit" /* 46 */
"implicit", /* 46 */
"same_model" /* 47 */
};
int count_opt_list = 47;
int count_opt_list = 48;
strcpy(file_name, "phreeqc.dmp");
/*
@ -113,7 +114,7 @@ read_transport(void)
}
else
old_cells = count_cells;
count_length = count_disp = count_punch = count_print = count_por = 0;
count_length = count_disp = count_punch = count_print = count_por = count_same_model = 0;
length = (LDBLE *)PHRQ_malloc(sizeof(LDBLE));
if (length == NULL)
@ -135,6 +136,10 @@ read_transport(void)
if (print_temp == NULL)
malloc_error();
same_model_temp = (int *)PHRQ_malloc(sizeof(int));
if (same_model_temp == NULL)
malloc_error();
count_length_alloc = count_disp_alloc = count_por_alloc = 1;
transport_start = 1;
/*
@ -706,8 +711,27 @@ read_transport(void)
//warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion.");
max_mixf = 1.0;
}
min_dif_LM = -30.0;
if (copy_token(token, &next_char, &l) != EMPTY)
{
/* minimal moles for diffusion */
if (sscanf(token, SCANFORMAT, &min_dif_LM) != 1)
{
input_error++;
error_string = sformatf(
"Expected the minimal log10(molality) for including a species in multicomponent diffusion,\n taking -30.0");
warning_msg(error_string);
break;
}
}
opt_save = OPTION_DEFAULT;
break;
case 47: /* same_model */
same_model_temp =
read_list_ints_range(&next_char, &count_same_model, FALSE,
same_model_temp);
opt_save = 47;
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
@ -767,6 +791,7 @@ read_transport(void)
cell_data[i].potV = 0;
cell_data[i].punch = FALSE;
cell_data[i].print = FALSE;
cell_data[i].same_model = FALSE;
}
all_cells = all_cells_now;
}
@ -961,6 +986,29 @@ read_transport(void)
else if (simul_tr == 1 || old_cells != count_cells)
for (i = 0; i < all_cells; i++)
cell_data[i].print = TRUE;
/*
* Fill in data for same_model
*/
if (count_same_model != 0)
{
for (i = 0; i < all_cells; i++)
cell_data[i].same_model = FALSE;
for (i = 0; i < count_same_model; i++)
{
if (same_model_temp[i] > all_cells - 1 || same_model_temp[i] < 0)
{
error_string = sformatf(
"Cell number for same_model is out of range, %d. Request ignored.",
same_model_temp[i]);
warning_msg(error_string);
}
else
cell_data[same_model_temp[i]].same_model = TRUE;
}
}
else if (simul_tr == 1 || old_cells != count_cells)
for (i = 0; i < all_cells; i++)
cell_data[i].same_model = FALSE;
//#define OLD_POROSITY
#if defined(OLD_POROSITY)
/*
@ -1043,7 +1091,7 @@ read_transport(void)
"Must enter diffusion coefficient (-diffc) when modeling thermal diffusion.");
error_msg(error_string, CONTINUE);
}
else if (heat_diffc > diffc)
else if (heat_diffc > diffc && !implicit)
{
error_string = sformatf(
"Thermal diffusion is calculated assuming exchange factor was for\n\t effective (non-thermal) diffusion coefficient = %e.",
@ -1053,7 +1101,7 @@ read_transport(void)
}
else
{
if (heat_diffc > diffc)
if (heat_diffc > diffc && !implicit)
{
input_error++;
error_string = sformatf(
@ -1077,6 +1125,7 @@ read_transport(void)
pors = (LDBLE *)free_check_null(pors);
punch_temp = (int *)free_check_null(punch_temp);
print_temp = (int *)free_check_null(print_temp);
same_model_temp = (int *)free_check_null(same_model_temp);
if (dump_in == TRUE)
{

1189
phreeqcpp/transport.cpp
View File

File diff suppressed because it is too large Load Diff