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

Merge commit 'aa34ac11538c27f39c44c34938411f1e00afce19'

Browse Source
This commit is contained in:
Scott R Charlton 2019-06-05 20:44:38 -06:00
commit 5a43b70d11
14 changed files with 1184 additions and 253 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
phreeqcpp/Makefile.old
View File

@ -971,8 +971,8 @@ tester:
# cd ../mytest; make clean; make -k -j 1 $(SPOOL) make.out $(SPOOL2); make diff $(SPOOL) diff.out $(SPOOL2)
cd ../mytest $(CONCAT) make clean $(CONCAT) make -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make diff $(SPOOL) diff.out $(SPOOL2)
cd ../examples $(CONCAT) make -f Makefile.old clean $(CONCAT) make -f Makefile.old -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make -f Makefile.old diff $(SPOOL) diff.out $(SPOOL2)
svn status -q ../mytest
svn status -q ../examples
# svn status -q ../mytest
# svn status -q ../examples
#ld-option
# Usage: ldflags += $(call ld-option, -Wl$(comma)--hash-style=sysv)

2
phreeqcpp/Phreeqc.cpp
View File

@ -711,6 +711,8 @@ void Phreeqc::init(void)
all_cells = 0;
multi_Dflag = FALSE;
interlayer_Dflag = FALSE;
implicit = FALSE;
max_mixf = 1.0;
default_Dw = 0;
correct_Dw = 0;
multi_Dpor = 0;

7
phreeqcpp/Phreeqc.h
View File

@ -1065,10 +1065,11 @@ public:
LDBLE viscosity(void);
LDBLE calc_vm_Cl(void);
int multi_D(LDBLE DDt, int mobile_cell, int stagnant);
int 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);
int fill_spec(int cell_no);
void define_ct_structures(void);
int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec);
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);
int diff_stag_surf(int mobile_cell);
@ -1426,6 +1427,8 @@ protected:
int old_cells, max_cells, all_cells;
int multi_Dflag; /* signals calc'n of multicomponent diffusion */
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 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 */

3
phreeqcpp/basicsubs.cpp
View File

@ -1470,7 +1470,8 @@ kinetics_moles(const char *kinetics_name)
error_string = sformatf( "No data for rate %s in KINETICS keyword.",
kinetics_name);
warning_msg(error_string);
//if (count_warnings >= 0) // appt debug cvode
// warning_msg(error_string);
return (0);
}
/* ---------------------------------------------------------------------- */

10
phreeqcpp/class_main.cpp
View File

@ -277,8 +277,8 @@ write_banner(void)
/* version */
#ifdef NPP
//len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8");
//len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1 AmpŠre");
len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1");
#else
len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
#endif
@ -302,7 +302,7 @@ write_banner(void)
/* date */
#ifdef NPP
len = sprintf(buffer, "%s", "January 17, 2019");
len = sprintf(buffer, "%s", "May 29, 2019");
#else
len = sprintf(buffer, "%s", "@VER_DATE@");
#endif
@ -492,12 +492,12 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
screen_msg(sformatf("Database file: %s\n\n", token));
strcpy(db_file, token);
#ifdef NPP
//output_msg(sformatf("Using PHREEQC: version 3.4.8 Ampère, compiled January 17, 2019\n"));
//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.4.8, compiled January 17, 2019\n"));
output_msg(sformatf("Using PHREEQC: version 3.5.1, compiled May 29, 2019\n"));
#endif
output_msg(sformatf("Database file: %s\n\n", token));
#ifdef NPP

4
phreeqcpp/common/PHRQ_io.cpp
View File

@ -32,7 +32,7 @@ PHRQ_io(void)
punch_on = true;
error_on = true;
dump_on = true;
echo_on = true; //**appt
echo_on = true;
screen_on = true;
echo_destination = ECHO_OUTPUT;
@ -800,7 +800,7 @@ get_line(void)
this->push_istream(next_stream);
std::ostringstream errstr;
errstr << "\n\tReading data from " << file_name <<" ...\n";
output_msg(errstr.str().c_str()); // **appt
output_msg(errstr.str().c_str());
}
continue;
}

4
phreeqcpp/global_structures.h
View File

@ -1052,12 +1052,12 @@ struct sol_D
struct J_ij
{
const char *name;
LDBLE tot1, tot2; /* species change in cells i and j */
LDBLE tot1, tot2, charge; /* species change in cells i and j */
};
struct M_S
{
const char *name;
LDBLE tot1, tot2;
LDBLE tot1, tot2, charge; /* master species transport in cells i and j */
};
// Pitzer definitions
typedef enum

57
phreeqcpp/kinetics.cpp
View File

@ -2478,7 +2478,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
* Rates and moles of each reaction are calculated in calc_kinetic_reaction
* Total number of moles in reaction is stored in kinetics[i].totals
*/
//int increase_tol = 0; // appt
int converge, m_iter;
int pr_all_save;
int nsaver;
@ -2529,9 +2529,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
converge =
set_and_run_wrapper(i, use_mix, FALSE, nsaver, step_fraction);
if (converge == MASS_BALANCE)
error_msg
("Negative concentration in system. Stopping calculation.",
STOP);
{
error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no);
error_msg(error_string, STOP);
}
run_reactions_iterations += iterations;
}
else
@ -2603,9 +2604,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
else
converge = set_and_run_wrapper(i, use_mix, FALSE, i, step_fraction);
if (converge == MASS_BALANCE)
error_msg
("Negative concentration in system. Stopping calculation.",
STOP);
{
error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no);
error_msg(error_string, STOP);
}
saver();
pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, i);
ss_assemblage_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, i);
@ -2675,6 +2677,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
cvode_mem = CVodeMalloc(n_reactions, f, 0.0, y, ADAMS, FUNCTIONAL, SV, &reltol, abstol, NULL, NULL, FALSE, iopt, ropt, machEnv);
iopt[MXSTEP] is maximum number of steps that CVODE tries.
*/
//iopt[SLDET] = TRUE; // appt
iopt[MXSTEP] = kinetics_ptr->Get_cvode_steps();
iopt[MAXORD] = kinetics_ptr->Get_cvode_order();
kinetics_cvode_mem =
@ -2704,30 +2707,49 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
*/
m_iter = 0;
sum_t = 0;
RESTART:
RESTART:
while (flag != SUCCESS)
{
sum_t += cvode_last_good_time;
if (state != TRANSPORT)
{
error_string = sformatf(
"CVode incomplete at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d, continuing...\n",
(int)iopt[NST], cell_no, (double)sum_t, m_iter + 1);
warning_msg(error_string);
error_string = sformatf("CV_ODE: Time: %8.2e s. Delta t: %8.2e s. Calls: %d.", (double)(sum_t), (double) cvode_last_good_time, m_iter);
status(0, error_string, true);
}
//if (state != TRANSPORT)
//{
// error_string = sformatf(
// "CVode incomplete at cvode_steps %d. Cell: %d. Time: %8.2e s. Cvode calls: %d, continuing...\n",
// (int)iopt[NST], cell_no, (double)sum_t, m_iter + 1);
// warning_msg(error_string);
//}
#ifdef DEBUG_KINETICS
if (m_iter > 5)
dump_kinetics_stderr(cell_no);
#endif
//if (m_iter > 0.5 * kinetics_ptr->Get_bad_step_max() &&
// (cvode_last_good_time < 1e-6 || cvode_last_good_time < 1e-6 * tout)) // appt
//{
// if (increase_tol < 3)
// {
// increase_tol += 1;
// for (size_t j = 0; j < kinetics_ptr->Get_kinetics_comps().size(); j++)
// {
// cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
// LDBLE tr = kinetics_comp_ptr->Get_tol() * 10.0;
// kinetics_comp_ptr->Set_tol(tr);
// tr += 0;
// }
// }
//}
cvode_last_good_time = 0;
if (++m_iter >= kinetics_ptr->Get_bad_step_max())
{
m_temp = (LDBLE *) free_check_null(m_temp);
m_original = (LDBLE *) free_check_null(m_original);
error_string = sformatf(
"Restart of integration at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d.",
(int)iopt[NST], cell_no, (double)sum_t, m_iter - 1);
"CVode is at maximum calls: %d. Cell: %d. Time: %8.2e s\nERROR: Please increase the maximum calls with -bad_step_max.",
m_iter, cell_no, (double)sum_t);
error_msg(error_string, STOP);
}
tout1 = tout - sum_t;
@ -2810,6 +2832,11 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction)
free_cvode();
use.Set_mix_in(use_save.Get_mix_in());
use.Set_mix_ptr(use_save.Get_mix_ptr());
error_string = sformatf("CV_ODE: Final Delta t: %8.2e s. Calls: %d. ", (double)cvode_last_good_time, m_iter);
status(0, error_string, true);
//status(0, NULL);
}
rate_sim_time = rate_sim_time_start + kin_time;

1
phreeqcpp/phqalloc.cpp
View File

@ -142,6 +142,7 @@ 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
if (p == NULL)
return NULL;

3
phreeqcpp/pitzer.cpp
View File

@ -2555,7 +2555,6 @@ gammas_pz(bool exch_a_f)
case 9: /* activity water */
break;
case 4: /* Exchange */
/*
* Find CEC
* z contains valence of cation for exchange species, alk contains cec
@ -2580,7 +2579,6 @@ gammas_pz(bool exch_a_f)
/*
* All other species
*/
/* modific 29 july 2005... */
if (s_x[i]->equiv != 0 && s_x[i]->alk > 0)
{
@ -2595,7 +2593,6 @@ gammas_pz(bool exch_a_f)
continue;
coef = s_x[i]->rxn_x->token[j].coef;
s_x[i]->lg += coef * s_x[i]->rxn_x->token[j].s->lg;
/*s_x[i]->dg += coef * s_x[i]->rxn_x->token[j].s->dg;*//* remove? */
}
}
if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles)

30
phreeqcpp/readtr.cpp
View File

@ -94,9 +94,10 @@ read_transport(void)
"porosities", /* 42 */
"porosity", /* 43 */
"fix_current", /* 44 */
"current" /* 45 */
"current", /* 45 */
"implicit" /* 46 */
};
int count_opt_list = 46;
int count_opt_list = 47;
strcpy(file_name, "phreeqc.dmp");
/*
@ -682,6 +683,31 @@ read_transport(void)
}
opt_save = OPTION_DEFAULT;
break;
case 46: /* implicit diffusion */
copy_token(token, &next_char, &l);
str_tolower(token);
if (strstr(token, "f") == token)
implicit = FALSE;
else if (strstr(token, "t") == token)
implicit = TRUE;
else
{
input_error++;
error_msg
("Expected flag for implicit diffusion calc`s: 'true' or 'false'.",
CONTINUE);
}
if (copy_token(token, &next_char, &l) == DIGIT)
{
sscanf(token, SCANFORMAT, &max_mixf);
}
else
{
//warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion.");
max_mixf = 1.0;
}
opt_save = OPTION_DEFAULT;
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;

2
phreeqcpp/tidy.cpp
View File

@ -1437,7 +1437,7 @@ tidy_inverse(void)
s_eminus->primary->in = TRUE; /* Include electrons */
if (master_alk_ptr)
{
master_alk_ptr->in = TRUE; /* Include alkalinity */
master_alk_ptr->in = TRUE; /* Include alkalinity */
}
else
{

1307
phreeqcpp/transport.cpp
View File

File diff suppressed because it is too large Load Diff

3
phreeqcpp/utilities.cpp
View File

@ -6,7 +6,6 @@
#include "Solution.h"
#include <time.h>
/* ---------------------------------------------------------------------- */
int Phreeqc::
add_elt_list(struct elt_list *elt_list_ptr, LDBLE coef)
@ -1516,7 +1515,7 @@ status(int count, const char *str, bool rk_string)
}
else
{
screen_string = sformatf("%-15s%-27s%1s%37s", sim_str, state_str, spin_str, stdstr.c_str());
screen_string = sformatf("%-15s%-27s%1s%45s", sim_str, state_str, spin_str, stdstr.c_str());
status_string = screen_string;
}
}