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

Merge commit '3fd8155d566798e121a1f54418f901cdd0db5e6c'

Browse Source
This commit is contained in:
Scott R Charlton 2019-01-28 21:20:14 -07:00
commit 43386a0bc4
3 changed files with 131 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
phreeqcpp/kinetics.cpp
View File

@ -1266,7 +1266,7 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
std::auto_ptr<cxxSSassemblage> ss_assemblage_save(NULL); std::auto_ptr<cxxSSassemblage> ss_assemblage_save(NULL);
std::auto_ptr<cxxKinetics> kinetics_save(NULL); std::auto_ptr<cxxKinetics> kinetics_save(NULL);
#endif #endif
int restart = 0;
small_pe_step = 5.; small_pe_step = 5.;
small_step = 10.; small_step = 10.;
@ -1312,14 +1312,16 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
{ {
diagonal_scale = TRUE; diagonal_scale = TRUE;
always_full_pitzer = FALSE; always_full_pitzer = FALSE;
max_try = 13; max_try = 15;
} }
else else
{ {
max_try = 14; max_try = 15;
} }
max_try = (max_tries < max_try) ? max_tries : max_try; max_try = (max_tries < max_try) ? max_tries : max_try;
/*max_try = 1; */ /*max_try = 1; */
restart:
for (j = 0; j < max_try; j++) for (j = 0; j < max_try; j++)
{ {
if (j == 1) if (j == 1)
@ -1468,6 +1470,26 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
error_string = sformatf( "Trying reduced tolerance %g ...\n", error_string = sformatf( "Trying reduced tolerance %g ...\n",
(double) ineq_tol); (double) ineq_tol);
warning_msg(error_string); warning_msg(error_string);
}
else if (j == 14 && use.Get_ss_assemblage_in())
{
//cxxStorageBin error_bin;
//Use2cxxStorageBin(error_bin);
//std::ostringstream error_input;
//error_bin.dump_raw(error_input, 0);
//cxxStorageBin reread;
//std::istringstream is(error_input.str());
//CParser cp(is);
//cp.set_echo_stream(CParser::EO_NONE);
//reread.read_raw(cp);
//cxxStorageBin2phreeqc(reread);
//error_string = sformatf("Trying restarting ...\n");
//warning_msg(error_string);
//if (restart < 2)
//{
// restart++;
// goto restart;
//}
} }
if (j > 0) if (j > 0)
{ {
@ -1487,6 +1509,29 @@ set_and_run_wrapper(int i, int use_mix, int use_kinetics, int nsaver,
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, kinetics_save->Get_n_user())); use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, kinetics_save->Get_n_user()));
} }
} }
if (j == 14)
{
cxxStorageBin error_bin;
Use2cxxStorageBin(error_bin);
std::ostringstream error_input;
error_bin.dump_raw(error_input, 0);
cxxStorageBin reread;
std::istringstream is(error_input.str());
CParser cp(is);
cp.set_echo_stream(CParser::EO_NONE);
reread.read_raw(cp);
cxxStorageBin2phreeqc(reread);
error_string = sformatf("Trying restarting ...\n");
warning_msg(error_string);
step_size = 1.0 + (small_step - 1.0)/((double) restart + 1.0);
pe_step_size = 1.0 + (small_pe_step - 1)/ ((double)restart + 1.0);
if (restart < 2)
{
restart++;
goto restart;
}
}
set_and_run_attempt = j; set_and_run_attempt = j;
converge = converge =

114
phreeqcpp/model.cpp
View File

@ -532,15 +532,15 @@ check_residuals(void)
{ {
if (x[i]->ss_in == FALSE) if (x[i]->ss_in == FALSE)
continue; continue;
if (x[i]->moles <= 1e0*MIN_TOTAL) //if (x[i]->moles <= 1e2*MIN_TOTAL)
continue; // continue;
if (residual[i] >= epsilon if (residual[i] >= epsilon
|| residual[i] <= -epsilon /* || stop_program == TRUE */ ) || residual[i] <= -epsilon /* || stop_program == TRUE */ )
{ {
error_string = sformatf( error_string = sformatf(
"%20s Total moles in solid solution has not converged. " "%20s Total moles in solid solution has not converged. "
"\tResidual: %e\n", x[i]->description, "\tResidual: %e %e\n", x[i]->description,
(double) residual[i]); (double) residual[i], x[i]->moles);
error_msg(error_string, CONTINUE); error_msg(error_string, CONTINUE);
} }
} }
@ -1696,7 +1696,7 @@ ineq(int in_kode)
*/ */
for (i = 0; i < count_unknowns; i++) for (i = 0; i < count_unknowns; i++)
{ {
if ((x[i]->type == PP || x[i]->type == SS_MOLES) if ((x[i]->type == PP || x[i]->type == SS_MOLES) && x[i]->phase->in == TRUE
&& pp_column_scale != 1.0) && pp_column_scale != 1.0)
{ {
for (j = 0; j < l_count_rows; j++) for (j = 0; j < l_count_rows; j++)
@ -2219,12 +2219,20 @@ mb_ss(void)
{ {
cxxSS *ss_ptr = ss_ptrs[i]; cxxSS *ss_ptr = ss_ptrs[i];
total_moles = 0; total_moles = 0;
bool ss_in = true;
for (size_t j = 0; j < ss_ptr->Get_ss_comps().size(); j++) for (size_t j = 0; j < ss_ptr->Get_ss_comps().size(); j++)
{ {
int l;
struct phase *phase_ptr = phase_bsearch(ss_ptr->Get_ss_comps()[j].Get_name().c_str(), &l, FALSE);
if (phase_ptr->in == FALSE)
{
continue;
}
cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[j]); cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[j]);
total_moles += comp_ptr->Get_moles(); total_moles += comp_ptr->Get_moles();
} }
if (total_moles > 1e-13) if (total_moles > 1.e10*MIN_TOTAL)
{ {
ss_ptr->Set_ss_in(true); ss_ptr->Set_ss_in(true);
} }
@ -2236,7 +2244,7 @@ mb_ss(void)
/* /*
* Calculate IAPc and IAPb * Calculate IAPc and IAPb
*/ */
if (phase0_ptr->rxn_x != NULL) if (phase0_ptr->in == TRUE && phase0_ptr->rxn_x != NULL)
{ {
log10_iap = 0; log10_iap = 0;
for (rxn_ptr = phase0_ptr->rxn_x->token + 1; for (rxn_ptr = phase0_ptr->rxn_x->token + 1;
@ -2250,7 +2258,7 @@ mb_ss(void)
{ {
iapc = 1e-99; iapc = 1e-99;
} }
if (phase1_ptr->rxn_x != NULL) if (phase1_ptr->in == TRUE && phase1_ptr->rxn_x != NULL)
{ {
log10_iap = 0; log10_iap = 0;
for (rxn_ptr = phase1_ptr->rxn_x->token + 1; for (rxn_ptr = phase1_ptr->rxn_x->token + 1;
@ -2340,9 +2348,12 @@ mb_ss(void)
{ {
if (x[i]->type != SS_MOLES) if (x[i]->type != SS_MOLES)
break; break;
//cxxSS *ss_ptr = use.Get_ss_assemblage_ptr()->Find(x[i]->ss_name);
cxxSS *ss_ptr = (cxxSS *) x[i]->ss_ptr; cxxSS *ss_ptr = (cxxSS *) x[i]->ss_ptr;
x[i]->ss_in = ss_ptr->Get_ss_in() ? TRUE : FALSE; x[i]->ss_in = FALSE;
if (x[i]->phase->in == TRUE && ss_ptr->Get_ss_in())
{
x[i]->ss_in = TRUE;
}
} }
return (OK); return (OK);
} }
@ -3213,7 +3224,7 @@ reset(void)
*/ */
for (i = 0; i < count_unknowns; i++) for (i = 0; i < count_unknowns; i++)
{ {
if (x[i]->type == PP || x[i]->type == SS_MOLES) if ((x[i]->type == PP || x[i]->type == SS_MOLES) && x[i]->phase->in == TRUE)
{ {
if (delta[i] < -1e8) if (delta[i] < -1e8)
@ -3232,21 +3243,21 @@ reset(void)
assert(comp_ptr); assert(comp_ptr);
if ((delta[i] < 0.0) if ((delta[i] < 0.0)
&& (-delta[i] > && (-delta[i] >
(comp_ptr->Get_initial_moles() - x[i]->moles))) (comp_ptr->Get_initial_moles() - x[i]->moles)))
{ {
if ((comp_ptr->Get_initial_moles() - x[i]->moles) != if ((comp_ptr->Get_initial_moles() - x[i]->moles) !=
0.0) 0.0)
{ {
f0 = fabs(delta[i] / f0 = fabs(delta[i] /
(comp_ptr->Get_initial_moles() - (comp_ptr->Get_initial_moles() -
x[i]->moles)); x[i]->moles));
if (f0 > factor) if (f0 > factor)
{ {
if (debug_model == TRUE) if (debug_model == TRUE)
{ {
output_msg(sformatf( output_msg(sformatf(
"%-10.10s, Precipitating too much dissolve_only mineral.\tDelta %e\tCurrent %e\tInitial %e\n", "%-10.10s, Precipitating too much dissolve_only mineral.\tDelta %e\tCurrent %e\tInitial %e\n",
x[i]->description, x[i]->description,
(double) delta[i], (double) delta[i],
(double) x[i]->moles, (double) x[i]->moles,
(double) comp_ptr->Get_initial_moles())); (double) comp_ptr->Get_initial_moles()));
@ -3259,9 +3270,9 @@ reset(void)
if (debug_model == TRUE) if (debug_model == TRUE)
{ {
output_msg(sformatf( output_msg(sformatf(
"%-10.10s, Precipitating dissolve_only mineral.\tDelta %e\n", "%-10.10s, Precipitating dissolve_only mineral.\tDelta %e\n",
x[i]->description, x[i]->description,
(double) delta[i])); (double) delta[i]));
} }
delta[i] = 0; delta[i] = 0;
} }
@ -3276,7 +3287,7 @@ reset(void)
if (debug_model == TRUE) if (debug_model == TRUE)
{ {
output_msg(sformatf( output_msg(sformatf(
"%-10.10s, Removing more than total mineral.\t%f\n", "%-10.10s, Removing more than total mineral.\t%f\n",
x[i]->description, (double) f0)); x[i]->description, (double) f0));
} }
factor = f0; factor = f0;
@ -3287,17 +3298,17 @@ reset(void)
if (debug_model == TRUE) if (debug_model == TRUE)
{ {
output_msg(sformatf( output_msg(sformatf(
"%-10.10s\tDelta: %e\tMass: %e " "%-10.10s\tDelta: %e\tMass: %e "
"Dissolving mineral with 0.0 mass.\n ", "Dissolving mineral with 0.0 mass.\n ",
x[i]->description, (double) delta[i], x[i]->description, (double) delta[i],
(double) x[i]->moles)); (double) x[i]->moles));
} }
delta[i] = 0.0; delta[i] = 0.0;
} }
else if (x[i]->ss_comp_name != NULL && delta[i] < -x[i]->phase->delta_max) else if (x[i]->ss_comp_name != NULL && delta[i] < -x[i]->phase->delta_max)
// Uses delta_max computed in step // Uses delta_max computed in step
// delta_max is the maximum amount of the mineral that could form based // delta_max is the maximum amount of the mineral that could form based
// on the limiting element in the system // on the limiting element in the system
{ {
f0 = -delta[i] / x[i]->phase->delta_max; f0 = -delta[i] / x[i]->phase->delta_max;
if (f0 > factor) if (f0 > factor)
@ -3305,24 +3316,38 @@ reset(void)
if (debug_model == TRUE) if (debug_model == TRUE)
{ {
output_msg(sformatf( output_msg(sformatf(
"%-10.10s, Precipitating too much mineral.\t%f\n", "%-10.10s, Precipitating too much mineral.\t%f\n",
x[i]->description, (double) f0)); x[i]->description, (double) f0));
} }
factor = f0; factor = f0;
} }
if (x[i]->moles > 0 && -delta[i] > 1e1*x[i]->moles) }
else if (x[i]->ss_comp_name != NULL && delta[i] < -pe_step_size*x[i]->moles)
{
f0 = (-delta[i]) / (pe_step_size*x[i]->moles);
if (f0 > factor)
{ {
f0 = (-delta[i])/(1e1*x[i]->moles); if (debug_model == TRUE)
if (f0 > factor)
{ {
if (debug_model == TRUE) output_msg(sformatf(
{ "%-10.10s, Precipitating too much mineral.\t%f\n",
output_msg(sformatf( x[i]->description, (double)f0));
"%-10.10s, Precipitating too much mineral.\t%f\n",
x[i]->description, (double)f0));
}
factor = f0;
} }
factor = f0;
}
}
else if (x[i]->ss_comp_name != NULL && delta[i] > x[i]->moles/ pe_step_size)
{
f0 = (delta[i]) / (x[i]->moles/ pe_step_size);
if (f0 > factor)
{
if (debug_model == TRUE)
{
output_msg(sformatf(
"%-10.10s, Precipitating too much mineral.\t%f\n",
x[i]->description, (double)f0));
}
factor = f0;
} }
} }
} }
@ -3908,7 +3933,7 @@ reset(void)
} }
last_patm_x = patm_x; last_patm_x = patm_x;
} }
else if (x[i]->type == SS_MOLES) else if (x[i]->type == SS_MOLES && x[i]->ss_in == TRUE)
{ {
/*if (fabs(delta[i]) > epsilon) converge=FALSE; */ /*if (fabs(delta[i]) > epsilon) converge=FALSE; */
@ -4266,9 +4291,10 @@ residuals(void)
} }
else if (x[i]->type == SS_MOLES) else if (x[i]->type == SS_MOLES)
{ {
if (x[i]->ss_in == FALSE)
continue;
residual[i] = x[i]->f * LOG_10; residual[i] = x[i]->f * LOG_10;
if (x[i]->moles <= 1e0*MIN_TOTAL) continue; if (fabs(residual[i]) > l_toler)
if (fabs(residual[i]) > l_toler && x[i]->ss_in == TRUE)
{ {
if (print_fail) if (print_fail)
output_msg(sformatf( output_msg(sformatf(
@ -5550,9 +5576,9 @@ numerical_jacobian(void)
molalities(TRUE); molalities(TRUE);
mb_sums(); mb_sums();
residuals(); residuals();
/* /*
* Clear array, note residuals are in array[i, count_unknowns+1] * Clear array, note residuals are in array[i, count_unknowns+1]
*/ */
for (i = 0; i < count_unknowns; i++) for (i = 0; i < count_unknowns; i++)
{ {
@ -5561,10 +5587,10 @@ numerical_jacobian(void)
for (i = 1; i < count_unknowns; i++) for (i = 1; i < count_unknowns; i++)
{ {
memcpy((void *) &(my_array[i * (count_unknowns + 1)]), memcpy((void *) &(my_array[i * (count_unknowns + 1)]),
(void *) &(my_array[0]), (size_t)count_unknowns * sizeof(LDBLE)); (void *) &(my_array[0]), (size_t) count_unknowns * sizeof(LDBLE));
} }
base = (LDBLE *)PHRQ_malloc((size_t)count_unknowns * sizeof(LDBLE)); base = (LDBLE *) PHRQ_malloc((size_t) count_unknowns * sizeof(LDBLE));
if (base == NULL) if (base == NULL)
malloc_error(); malloc_error();
for (i = 0; i < count_unknowns; i++) for (i = 0; i < count_unknowns; i++)

14
phreeqcpp/structures.cpp
View File

@ -3746,11 +3746,23 @@ Use2cxxStorageBin(cxxStorageBin & sb)
sb.Get_system().Set_io(sb.Get_io()); sb.Get_system().Set_io(sb.Get_io());
if (use.Get_mix_in()) if (use.Get_mix_in())
{ {
cxxMix *entity = Utilities::Rxn_find(Rxn_mix_map, use.Get_n_mix_user()); cxxMix *entity = use.Get_mix_ptr();
if (entity != NULL) if (entity != NULL)
{ {
sb.Set_Mix(use.Get_n_mix_user(), entity); sb.Set_Mix(use.Get_n_mix_user(), entity);
} }
// put mix solutions in sb
cxxMix * mix_ptr = use.Get_mix_ptr();
std::map<int, LDBLE>::const_iterator cit;
for (cit = mix_ptr->Get_mixComps().begin(); cit != mix_ptr->Get_mixComps().end(); cit++)
{
cxxSolution *entity = Utilities::Rxn_find(Rxn_solution_map, cit->first);
if (entity != NULL)
{
sb.Set_Solution(cit->first, entity);
}
}
} }
else if (use.Get_solution_in()) else if (use.Get_solution_in())
{ {