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

initial Tony changes

Browse Source
This commit is contained in:
David Parkhurst 2019-10-22 16:36:15 -06:00
parent 009aec74d1
commit 8089c10273
15 changed files with 1030 additions and 363 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
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
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 */

4
SurfaceComp.cxx
View File

@ -414,7 +414,7 @@ cxxSurfaceComp::add(const cxxSurfaceComp & addee, LDBLE extensive)
{
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 +430,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
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
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
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

5
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

8
integrate.cpp
View File

@ -1069,6 +1069,14 @@ 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());

4
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);
}

4
mainsubs.cpp
View File

@ -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
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
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
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
*/

63
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)
{

1087
transport.cpp
View File

File diff suppressed because it is too large Load Diff