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

Merge commit '63108ac76d6fa7c378c21690b6342333c4cb446c'

Browse Source
This commit is contained in:
Scott R Charlton 2019-10-30 11:10:37 -06:00
commit 545e81ba3a
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->la = f1 * this->la + f2 * addee.la;
this->charge_balance += addee.charge_balance * extensive; 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; std::ostringstream oss;
oss << oss <<
@ -374,7 +375,8 @@ cxxExchComp::add(const cxxExchComp & addee, LDBLE extensive)
this->phase_proportion = this->phase_proportion =
this->phase_proportion * f1 + addee.phase_proportion * f2; 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; std::ostringstream oss;
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(); this->exchange_comps.clear();
for (int n = 0; n < count; n++) for (int n = 0; n < count; n++)
{ {
cxxExchComp ec; cxxExchComp ec(this->io);
ec.Deserialize(dictionary, ints, doubles, ii, dd); ec.Deserialize(dictionary, ints, doubles, ii, dd);
this->exchange_comps.push_back(ec); 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(); this->pp_assemblage_comps.clear();
for (int n = 0; n < count; n++) for (int n = 0; n < count; n++)
{ {
cxxPPassemblageComp ppc; cxxPPassemblageComp ppc(this->io);
ppc.Deserialize(dictionary, ints, doubles, ii, dd); ppc.Deserialize(dictionary, ints, doubles, ii, dd);
std::string str(ppc.Get_name()); std::string str(ppc.Get_name());
this->pp_assemblage_comps[str] = ppc; this->pp_assemblage_comps[str] = ppc;

2
phreeqcpp/Phreeqc.cpp
View File

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

9
phreeqcpp/Phreeqc.h
View File

@ -1066,9 +1066,11 @@ public:
LDBLE calc_vm_Cl(void); LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant); int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, 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); void diffuse_implicit(LDBLE DDt, int stagnant);
int fill_spec(int cell_no); int fill_spec(int cell_no, int ref_cell);
void define_ct_structures(void); 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); 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); static int sort_species_name(const void *ptr1, const void *ptr2);
int disp_surf(LDBLE stagkin_time); int disp_surf(LDBLE stagkin_time);
@ -1429,6 +1431,7 @@ protected:
int interlayer_Dflag; /* multicomponent diffusion and diffusion through interlayer porosity */ int interlayer_Dflag; /* multicomponent diffusion and diffusion through interlayer porosity */
int implicit; /* implicit calculation of diffusion */ int implicit; /* implicit calculation of diffusion */
LDBLE max_mixf; /* the maximum value of the implicit mixfactor = De * Dt / (Dx^2) */ 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 */ LDBLE default_Dw; /* default species diffusion coefficient in water at 25oC, m2/s */
int correct_Dw; /* if true, Dw is adapted in calc_SC */ int correct_Dw; /* if true, Dw is adapted in calc_SC */
LDBLE multi_Dpor; /* uniform porosity of free porewater in solid medium */ 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(); this->ss_comps.clear();
for (int i = 0; i < count; i++) for (int i = 0; i < count; i++)
{ {
cxxSScomp ssc; cxxSScomp ssc(this->io);
ssc.Deserialize(dictionary, ints, doubles, ii, dd); ssc.Deserialize(dictionary, ints, doubles, ii, dd);
this->ss_comps.push_back(ssc); 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(); this->SSs.clear();
for (int n = 0; n < count; n++) for (int n = 0; n < count; n++)
{ {
cxxSS ssc; cxxSS ssc(this->io);
ssc.Deserialize(dictionary, ints, doubles, ii, dd); ssc.Deserialize(dictionary, ints, doubles, ii, dd);
std::string str(ssc.Get_name()); std::string str(ssc.Get_name());
this->SSs[str] = ssc; this->SSs[str] = ssc;

6
phreeqcpp/Serializer.cxx
View File

@ -133,7 +133,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break; break;
case PT_GASPHASE: case PT_GASPHASE:
{ {
cxxGasPhase entity; cxxGasPhase entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd); entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user(); int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_gas_phase_map()[n_user] = entity; phreeqc_ref.Get_Rxn_gas_phase_map()[n_user] = entity;
@ -141,7 +141,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break; break;
case PT_KINETICS: case PT_KINETICS:
{ {
cxxKinetics entity; cxxKinetics entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd); entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user(); int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_kinetics_map()[n_user] = entity; phreeqc_ref.Get_Rxn_kinetics_map()[n_user] = entity;
@ -158,7 +158,7 @@ Serializer::Deserialize(Phreeqc &phreeqc_ref, Dictionary &dictionary, std::vecto
break; break;
case PT_SSASSEMBLAGE: case PT_SSASSEMBLAGE:
{ {
cxxSSassemblage entity; cxxSSassemblage entity(phreeqc_ref.Get_phrq_io());
entity.Deserialize(dictionary, ints, doubles, ii, dd); entity.Deserialize(dictionary, ints, doubles, ii, dd);
int n_user = entity.Get_n_user(); int n_user = entity.Get_n_user();
phreeqc_ref.Get_Rxn_ss_assemblage_map()[n_user] = entity; 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; break;
case Keywords::KEY_SOLID_SOLUTIONS_RAW: case Keywords::KEY_SOLID_SOLUTIONS_RAW:
{ {
cxxSSassemblage entity; cxxSSassemblage entity(this->Get_io());
entity.read_raw(parser); entity.read_raw(parser);
SSassemblages[entity.Get_n_user()] = entity; SSassemblages[entity.Get_n_user()] = entity;
} }
@ -1195,7 +1195,7 @@ cxxStorageBin::read_raw_keyword(CParser & parser)
case Keywords::KEY_SOLID_SOLUTIONS_RAW: case Keywords::KEY_SOLID_SOLUTIONS_RAW:
{ {
cxxSSassemblage entity; cxxSSassemblage entity(this->Get_io());
entity.read_raw(parser); entity.read_raw(parser);
SSassemblages[entity.Get_n_user()] = entity; SSassemblages[entity.Get_n_user()] = entity;
entity_number = entity.Get_n_user(); 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(); this->surface_comps.clear();
for (int n = 0; n < count; n++) for (int n = 0; n < count; n++)
{ {
cxxSurfaceComp sc; cxxSurfaceComp sc(this->io);
sc.Deserialize(dictionary, ints, doubles, ii, dd); sc.Deserialize(dictionary, ints, doubles, ii, dd);
this->surface_comps.push_back(sc); this->surface_comps.push_back(sc);
} }
@ -822,7 +822,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints,
this->surface_charges.clear(); this->surface_charges.clear();
for (int n = 0; n < count; n++) for (int n = 0; n < count; n++)
{ {
cxxSurfaceCharge sc; cxxSurfaceCharge sc(this->io);
sc.Deserialize(dictionary, ints, doubles, ii, dd); sc.Deserialize(dictionary, ints, doubles, ii, dd);
this->surface_charges.push_back(sc); 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->totals.add_extensive(addee.totals, extensive);
this->la = f1 * this->la + f2 * addee.la; this->la = f1 * this->la + f2 * addee.la;
this->charge_balance += addee.charge_balance * extensive; 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; std::ostringstream oss;
oss << oss <<
@ -410,11 +411,12 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
this->phase_proportion * f1 + addee.phase_proportion * f2; 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; std::ostringstream oss;
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; << this->formula;
error_msg(oss.str().c_str(), CONTINUE); error_msg(oss.str().c_str(), CONTINUE);
return; return;
@ -430,7 +432,7 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
{ {
std::ostringstream oss; std::ostringstream oss;
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; << this->formula;
error_msg(oss.str().c_str(), CONTINUE); error_msg(oss.str().c_str(), CONTINUE);
return; return;

10
phreeqcpp/class_main.cpp
View File

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

1
phreeqcpp/common/PHRQ_io.cpp
View File

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

81
phreeqcpp/cvode.cpp
View File

@ -72,6 +72,7 @@
#include "sundialsmath.h" #include "sundialsmath.h"
#include "Phreeqc.h" #include "Phreeqc.h"
#if !defined(WIN32_MEMORY_DEBUG) #if !defined(WIN32_MEMORY_DEBUG)
#define malloc MACHENV_MALLOC PHRQ_malloc #define malloc MACHENV_MALLOC PHRQ_malloc
#endif #endif
@ -1974,7 +1975,6 @@ CVStep(CVodeMem cv_mem)
} }
#endif #endif
CVSet(cv_mem); CVSet(cv_mem);
nflag = CVnls(cv_mem, nflag); nflag = CVnls(cv_mem, nflag);
if (CVMEM cvode_error == TRUE || predict_fail) if (CVMEM cvode_error == TRUE || predict_fail)
{ {
@ -2785,7 +2785,6 @@ CVnlsNewton(CVodeMem cv_mem, int nflag)
loop loop
{ {
f(N, tn, zn[0], ftemp, f_data); f(N, tn, zn[0], ftemp, f_data);
nfe++; nfe++;
@ -3408,7 +3407,7 @@ CVHandleFailure(CVodeMem cv_mem, int kflag)
size is reset accordingly. size is reset accordingly.
*****************************************************************/ *****************************************************************/
#ifdef ORIGINAL
void void
CVBDFStab(CVodeMem cv_mem) CVBDFStab(CVodeMem cv_mem)
{ {
@ -3469,6 +3468,82 @@ CVBDFStab(CVodeMem cv_mem)
nscon = 0; 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 ************************************ /********************* CVsldet ************************************
This routine detects stability limitation using stored scaled 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(); this->kinetics_comps.clear();
for (int i = 0; i < n; i++) for (int i = 0; i < n; i++)
{ {
cxxKineticsComp kc; cxxKineticsComp kc(this->io);
kc.Deserialize(dictionary, ints, doubles, ii, dd); kc.Deserialize(dictionary, ints, doubles, ii, dd);
this->kinetics_comps.push_back(kc); this->kinetics_comps.push_back(kc);
} }

5
phreeqcpp/global_structures.h
View File

@ -547,6 +547,7 @@ struct cell_data
LDBLE potV; /* potential (V) */ LDBLE potV; /* potential (V) */
int punch; int punch;
int print; int print;
int same_model;
}; };
/*---------------------------------------------------------------------- /*----------------------------------------------------------------------
@ -1052,12 +1053,12 @@ struct sol_D
struct J_ij struct J_ij
{ {
const char *name; 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 struct M_S
{ {
const char *name; 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 // Pitzer definitions
typedef enum typedef enum

12
phreeqcpp/integrate.cpp
View File

@ -1069,9 +1069,17 @@ calc_psi_avg(cxxSurfaceCharge *charge_ptr, LDBLE surf_chrg_eq)
l_iter++; l_iter++;
if (l_iter > 50) 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( error_string = sformatf(
"\nToo many iterations in subroutine calc_psi_avg; surface charge = %12.4e; surface water = %12.4e.\n", "\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()); (double)surf_chrg_eq, (double)charge_ptr->Get_mass_water());
error_msg(error_string, STOP); 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()) if (step_bad > kinetics_ptr->Get_bad_step_max())
{ {
error_string = sformatf( error_string = sformatf(
"Bad RK steps > %d. Please decrease (time)step or increase -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()); kinetics_ptr->Get_bad_step_max(), cell_no);
error_msg(error_string, STOP); error_msg(error_string, STOP);
} }
@ -1511,11 +1511,11 @@ restart:
} }
if (j == 14) if (j == 14)
{ {
cxxStorageBin error_bin; cxxStorageBin error_bin(this->Get_phrq_io());
Use2cxxStorageBin(error_bin); Use2cxxStorageBin(error_bin);
std::ostringstream error_input; std::ostringstream error_input;
error_bin.dump_raw(error_input, 0); error_bin.dump_raw(error_input, 0);
cxxStorageBin reread; cxxStorageBin reread(this->Get_phrq_io());
std::istringstream is(error_input.str()); std::istringstream is(error_input.str());
CParser cp(is); CParser cp(is);
cp.set_echo_stream(CParser::EO_NONE); cp.set_echo_stream(CParser::EO_NONE);
@ -1574,7 +1574,7 @@ restart:
* write to error.inp what failed to converge. * write to error.inp what failed to converge.
*/ */
std::ofstream error_input("error.inp"); std::ofstream error_input("error.inp");
cxxStorageBin error_bin; cxxStorageBin error_bin(this->Get_phrq_io());
Use2cxxStorageBin(error_bin); Use2cxxStorageBin(error_bin);
error_bin.dump_raw(error_input, 0); error_bin.dump_raw(error_input, 0);
error_input.close(); 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 * Save ss_assemblage composition into structure ss_assemblage with user
* number n_user. * number n_user.
*/ */
cxxSSassemblage temp_ss_assemblage; cxxSSassemblage temp_ss_assemblage(this->phrq_io);
if (use.Get_ss_assemblage_ptr() == NULL) if (use.Get_ss_assemblage_ptr() == NULL)
return (OK); return (OK);
@ -1719,6 +1719,8 @@ xsurface_save(int n_user)
if (x[i]->type == SURFACE) if (x[i]->type == SURFACE)
{ {
cxxSurfaceComp *comp_ptr = temp_surface.Find_comp(x[i]->surface_comp); 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); assert(comp_ptr);
comp_ptr->Set_la(x[i]->master[0]->s->la); comp_ptr->Set_la(x[i]->master[0]->s->la);
comp_ptr->Set_moles(0.); 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)) 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); 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); assert(charge_ptr);
charge_ptr->Set_charge_balance(x[i]->f); charge_ptr->Set_charge_balance(x[i]->f);
charge_ptr->Set_la_psi(x[i]->master[0]->s->la); 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)); assert((s_pTail == NULL) || (s_pTail->pNext == NULL));
p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num); //p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num);
//p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt
if (p == NULL) if (p == NULL)
return NULL; return NULL;

2
phreeqcpp/prep.cpp
View File

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

5
phreeqcpp/print.cpp
View File

@ -1353,8 +1353,9 @@ print_pp_assemblage(void)
{ {
if (x[j]->type != PP) if (x[j]->type != PP)
continue; continue;
//cxxPPassemblageComp * comp_ptr = pp_assemblage_ptr->Find(x[j]->pp_assemblage_comp_name); cxxPPassemblage * pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, use.Get_n_pp_assemblage_user());
cxxPPassemblageComp * comp_ptr = (cxxPPassemblageComp * ) x[j]->pp_assemblage_comp_ptr; 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 * Print saturation index
*/ */

13
phreeqcpp/read.cpp
View File

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

63
phreeqcpp/readtr.cpp
View File

@ -37,13 +37,13 @@ read_transport(void)
*/ */
char *ptr; char *ptr;
int i, j, l; 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; int count_length_alloc, count_disp_alloc, count_por_alloc;
char token[MAX_LENGTH]; char token[MAX_LENGTH];
char *description; char *description;
int n_user, n_user_end; int n_user, n_user_end;
LDBLE *length, *disp, *pors; LDBLE *length, *disp, *pors;
int *punch_temp, *print_temp; int *punch_temp, *print_temp, *same_model_temp;
int return_value, opt, opt_save; int return_value, opt, opt_save;
char *next_char, *next_char_save; char *next_char, *next_char_save;
char file_name[MAX_LENGTH]; char file_name[MAX_LENGTH];
@ -95,9 +95,10 @@ read_transport(void)
"porosity", /* 43 */ "porosity", /* 43 */
"fix_current", /* 44 */ "fix_current", /* 44 */
"current", /* 45 */ "current", /* 45 */
"implicit" /* 46 */ "implicit", /* 46 */
"same_model" /* 47 */
}; };
int count_opt_list = 47; int count_opt_list = 48;
strcpy(file_name, "phreeqc.dmp"); strcpy(file_name, "phreeqc.dmp");
/* /*
@ -113,7 +114,7 @@ read_transport(void)
} }
else else
old_cells = count_cells; 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)); length = (LDBLE *)PHRQ_malloc(sizeof(LDBLE));
if (length == NULL) if (length == NULL)
@ -135,6 +136,10 @@ read_transport(void)
if (print_temp == NULL) if (print_temp == NULL)
malloc_error(); 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; count_length_alloc = count_disp_alloc = count_por_alloc = 1;
transport_start = 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."); //warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion.");
max_mixf = 1.0; 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; opt_save = OPTION_DEFAULT;
break; 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) if (return_value == EOF || return_value == KEYWORD)
break; break;
@ -767,6 +791,7 @@ read_transport(void)
cell_data[i].potV = 0; cell_data[i].potV = 0;
cell_data[i].punch = FALSE; cell_data[i].punch = FALSE;
cell_data[i].print = FALSE; cell_data[i].print = FALSE;
cell_data[i].same_model = FALSE;
} }
all_cells = all_cells_now; all_cells = all_cells_now;
} }
@ -961,6 +986,29 @@ read_transport(void)
else if (simul_tr == 1 || old_cells != count_cells) else if (simul_tr == 1 || old_cells != count_cells)
for (i = 0; i < all_cells; i++) for (i = 0; i < all_cells; i++)
cell_data[i].print = TRUE; 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 //#define OLD_POROSITY
#if defined(OLD_POROSITY) #if defined(OLD_POROSITY)
/* /*
@ -1043,7 +1091,7 @@ read_transport(void)
"Must enter diffusion coefficient (-diffc) when modeling thermal diffusion."); "Must enter diffusion coefficient (-diffc) when modeling thermal diffusion.");
error_msg(error_string, CONTINUE); error_msg(error_string, CONTINUE);
} }
else if (heat_diffc > diffc) else if (heat_diffc > diffc && !implicit)
{ {
error_string = sformatf( error_string = sformatf(
"Thermal diffusion is calculated assuming exchange factor was for\n\t effective (non-thermal) diffusion coefficient = %e.", "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 else
{ {
if (heat_diffc > diffc) if (heat_diffc > diffc && !implicit)
{ {
input_error++; input_error++;
error_string = sformatf( error_string = sformatf(
@ -1077,6 +1125,7 @@ read_transport(void)
pors = (LDBLE *)free_check_null(pors); pors = (LDBLE *)free_check_null(pors);
punch_temp = (int *)free_check_null(punch_temp); punch_temp = (int *)free_check_null(punch_temp);
print_temp = (int *)free_check_null(print_temp); print_temp = (int *)free_check_null(print_temp);
same_model_temp = (int *)free_check_null(same_model_temp);
if (dump_in == TRUE) if (dump_in == TRUE)
{ {

1189
phreeqcpp/transport.cpp
View File

File diff suppressed because it is too large Load Diff