mirror of
https://git.gfz-potsdam.de/naaice/iphreeqc.git
synced 2025-12-15 16:18:22 +01:00
refactor(initialization): relocate initial setup methods
This commit is contained in:
Max Lübke
parent
7806aad339
commit
643f53890e
@ -338,597 +338,6 @@ int Phreeqc::set_use(void)
|
||||
*/
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
initial_solutions(int print)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/*
|
||||
* Go through list of solutions, make initial solution calculations
|
||||
* for any marked "new".
|
||||
*/
|
||||
int converge, converge1;
|
||||
int last, n_user, print1;
|
||||
char token[2 * MAX_LENGTH];
|
||||
|
||||
state = INITIAL_SOLUTION;
|
||||
set_use();
|
||||
print1 = TRUE;
|
||||
dl_type_x = cxxSurface::NO_DL;
|
||||
//std::map<int, cxxSolution>::iterator it = Rxn_solution_map.begin();
|
||||
//for ( ; it != Rxn_solution_map.end(); it++)
|
||||
//{
|
||||
//for (size_t nn = 0; nn < Rxn_new_solution.size(); nn++)
|
||||
for (std::set<int>::const_iterator nit = Rxn_new_solution.begin(); nit != Rxn_new_solution.end(); nit++)
|
||||
{
|
||||
std::map<int, cxxSolution>::iterator it = Rxn_solution_map.find(*nit);
|
||||
if (it == Rxn_solution_map.end())
|
||||
{
|
||||
assert(false);
|
||||
}
|
||||
cxxSolution &solution_ref = it->second;
|
||||
initial_solution_isotopes = FALSE;
|
||||
if (solution_ref.Get_new_def())
|
||||
{
|
||||
if (print1 == TRUE && print == TRUE)
|
||||
{
|
||||
dup_print("Beginning of initial solution calculations.",
|
||||
TRUE);
|
||||
print1 = FALSE;
|
||||
}
|
||||
if (print == TRUE)
|
||||
{
|
||||
snprintf(token, sizeof(token), "Initial solution %d.\t%.350s",
|
||||
solution_ref.Get_n_user(), solution_ref.Get_description().c_str());
|
||||
dup_print(token, FALSE);
|
||||
}
|
||||
use.Set_solution_ptr(&solution_ref);
|
||||
LDBLE d0 = solution_ref.Get_density();
|
||||
//LDBLE d1 = 0;
|
||||
bool diag = (diagonal_scale == TRUE) ? true : false;
|
||||
int count_iterations = 0;
|
||||
std::string input_units = solution_ref.Get_initial_data()->Get_units();
|
||||
cxxISolution *initial_data_ptr = solution_ref.Get_initial_data();
|
||||
density_iterations = 0;
|
||||
for (;;)
|
||||
{
|
||||
prep();
|
||||
k_temp(solution_ref.Get_tc(), solution_ref.Get_patm());
|
||||
set(TRUE);
|
||||
always_full_pitzer = FALSE;
|
||||
|
||||
diagonal_scale = TRUE;
|
||||
converge = model();
|
||||
if (converge == ERROR /*&& diagonal_scale == FALSE*/)
|
||||
{
|
||||
diagonal_scale = TRUE;
|
||||
always_full_pitzer = TRUE;
|
||||
set(TRUE);
|
||||
converge = model();
|
||||
}
|
||||
calc_dens();
|
||||
kgw_kgs = mass_water_aq_x / solution_mass_x;
|
||||
density_iterations++;
|
||||
if (solution_ref.Get_initial_data()->Get_calc_density())
|
||||
{
|
||||
solution_ref.Set_density(calc_dens());
|
||||
if (!equal(d0, solution_ref.Get_density(), 1e-8))
|
||||
{
|
||||
initial_data_ptr->Set_units(input_units);
|
||||
d0 = solution_ref.Get_density();
|
||||
if (count_iterations++ < 20)
|
||||
{
|
||||
diag = (diagonal_scale == TRUE) ? true : false;
|
||||
continue;
|
||||
}
|
||||
else
|
||||
{
|
||||
error_msg(sformatf("%s %d.", "Density calculation failed for initial solution ", solution_ref.Get_n_user()),
|
||||
STOP);
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
diagonal_scale = (diag) ? TRUE : FALSE;
|
||||
converge1 = check_residuals();
|
||||
sum_species();
|
||||
viscos = viscosity(NULL);
|
||||
use.Get_solution_ptr()->Set_viscosity(viscos);
|
||||
use.Get_solution_ptr()->Set_viscos_0(viscos_0);
|
||||
if (use.Get_surface_ptr() != NULL && dl_type_x != cxxSurface::NO_DL)
|
||||
use.Get_surface_ptr()->Set_DDL_viscosity(viscosity(use.Get_surface_ptr()));
|
||||
add_isotopes(solution_ref);
|
||||
punch_all();
|
||||
print_all();
|
||||
density_iterations = 0;
|
||||
/* free_model_allocs(); */
|
||||
// remove pr_in
|
||||
for (size_t i = 0; i < count_unknowns; i++)
|
||||
{
|
||||
if (x[i]->type == SOLUTION_PHASE_BOUNDARY)
|
||||
x[i]->phase->pr_in = false;
|
||||
}
|
||||
|
||||
if (converge == ERROR || converge1 == ERROR)
|
||||
{
|
||||
error_msg(sformatf("%s %d.", "Model failed to converge for initial solution ", solution_ref.Get_n_user()),
|
||||
STOP);
|
||||
}
|
||||
n_user = solution_ref.Get_n_user();
|
||||
last = solution_ref.Get_n_user_end();
|
||||
/* copy isotope data */
|
||||
if (solution_ref.Get_isotopes().size() > 0)
|
||||
{
|
||||
isotopes_x = solution_ref.Get_isotopes();
|
||||
}
|
||||
else
|
||||
{
|
||||
isotopes_x.clear();
|
||||
}
|
||||
xsolution_save(n_user);
|
||||
Utilities::Rxn_copies(Rxn_solution_map, n_user, last);
|
||||
}
|
||||
}
|
||||
initial_solution_isotopes = FALSE;
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
initial_exchangers(int print)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/*
|
||||
* Go through list of exchangers, make initial calculations
|
||||
* for any marked "new" that are defined to be in equilibrium with a
|
||||
* solution.
|
||||
*/
|
||||
int i, converge, converge1;
|
||||
int last, n_user, print1;
|
||||
char token[2 * MAX_LENGTH];
|
||||
|
||||
state = INITIAL_EXCHANGE;
|
||||
set_use();
|
||||
print1 = TRUE;
|
||||
dl_type_x = cxxSurface::NO_DL;
|
||||
//std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.begin();
|
||||
//for ( ; it != Rxn_exchange_map.end(); it++)
|
||||
//{
|
||||
//for (size_t nn = 0; nn < Rxn_new_exchange.size(); nn++)
|
||||
//{
|
||||
//std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.find(Rxn_new_exchange[nn]);
|
||||
for (std::set<int>::const_iterator nit = Rxn_new_exchange.begin(); nit != Rxn_new_exchange.end(); nit++)
|
||||
{
|
||||
std::map<int, cxxExchange>::iterator it = Rxn_exchange_map.find(*nit);
|
||||
if (it == Rxn_exchange_map.end())
|
||||
{
|
||||
assert(false);
|
||||
}
|
||||
if (!it->second.Get_new_def())
|
||||
continue;
|
||||
cxxExchange *exchange_ptr = &(it->second);
|
||||
n_user = exchange_ptr->Get_n_user();
|
||||
last = exchange_ptr->Get_n_user_end();
|
||||
exchange_ptr->Set_n_user_end(n_user);
|
||||
exchange_ptr->Set_new_def(false);
|
||||
if (exchange_ptr->Get_solution_equilibria())
|
||||
{
|
||||
if (print1 == TRUE && print == TRUE)
|
||||
{
|
||||
dup_print("Beginning of initial exchange"
|
||||
"-composition calculations.", TRUE);
|
||||
print1 = FALSE;
|
||||
}
|
||||
if (print == TRUE)
|
||||
{
|
||||
snprintf(token, sizeof(token), "Exchange %d.\t%.350s",
|
||||
exchange_ptr->Get_n_user(), exchange_ptr->Get_description().c_str());
|
||||
dup_print(token, FALSE);
|
||||
}
|
||||
use.Set_exchange_ptr(exchange_ptr);
|
||||
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, exchange_ptr->Get_n_solution()));
|
||||
if (use.Get_solution_ptr() == NULL)
|
||||
{
|
||||
error_msg
|
||||
("Solution not found for initial exchange calculation",
|
||||
STOP);
|
||||
}
|
||||
|
||||
prep();
|
||||
k_temp(use.Get_solution_ptr()->Get_tc(), use.Get_solution_ptr()->Get_patm());
|
||||
set(TRUE);
|
||||
converge = model();
|
||||
converge1 = check_residuals();
|
||||
sum_species();
|
||||
species_list_sort();
|
||||
print_exchange();
|
||||
if (pr.user_print)
|
||||
print_user_print();
|
||||
xexchange_save(n_user);
|
||||
punch_all();
|
||||
/* free_model_allocs(); */
|
||||
if (converge == ERROR || converge1 == ERROR)
|
||||
{
|
||||
error_msg
|
||||
("Model failed to converge for initial exchange calculation.",
|
||||
STOP);
|
||||
}
|
||||
}
|
||||
for (i = n_user + 1; i <= last; i++)
|
||||
{
|
||||
Utilities::Rxn_copy(Rxn_exchange_map, n_user, i);
|
||||
}
|
||||
}
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
initial_gas_phases(int print)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/*
|
||||
* Go through list of gas_phases, make initial calculations
|
||||
* for any marked "new" that are defined to be in equilibrium with a
|
||||
* solution.
|
||||
*/
|
||||
int converge, converge1;
|
||||
int last, n_user, print1;
|
||||
char token[2 * MAX_LENGTH];
|
||||
class phase *phase_ptr;
|
||||
class rxn_token *rxn_ptr;
|
||||
LDBLE lp;
|
||||
bool PR = false;
|
||||
|
||||
state = INITIAL_GAS_PHASE;
|
||||
set_use();
|
||||
print1 = TRUE;
|
||||
dl_type_x = cxxSurface::NO_DL;
|
||||
//std::map<int, cxxGasPhase>::iterator it = Rxn_gas_phase_map.begin();
|
||||
//for ( ; it != Rxn_gas_phase_map.end(); it++)
|
||||
//{
|
||||
//for (size_t nn = 0; nn < Rxn_new_gas_phase.size(); nn++)
|
||||
//{
|
||||
//std::map<int, cxxGasPhase>::iterator it = Rxn_gas_phase_map.find(Rxn_new_gas_phase[nn]);
|
||||
for (std::set<int>::const_iterator nit = Rxn_new_gas_phase.begin(); nit != Rxn_new_gas_phase.end(); nit++)
|
||||
{
|
||||
std::map<int, cxxGasPhase>::iterator it = Rxn_gas_phase_map.find(*nit);
|
||||
if (it == Rxn_gas_phase_map.end())
|
||||
{
|
||||
assert(false);
|
||||
}
|
||||
cxxGasPhase *gas_phase_ptr = &it->second;
|
||||
if (!gas_phase_ptr->Get_new_def())
|
||||
continue;
|
||||
n_user = gas_phase_ptr->Get_n_user();
|
||||
last = gas_phase_ptr->Get_n_user_end();
|
||||
gas_phase_ptr->Set_n_user_end(n_user);
|
||||
gas_phase_ptr->Set_new_def(false);
|
||||
if (gas_phase_ptr->Get_solution_equilibria())
|
||||
{
|
||||
if (print1 == TRUE && print == TRUE)
|
||||
{
|
||||
dup_print("Beginning of initial gas_phase"
|
||||
"-composition calculations.", TRUE);
|
||||
print1 = FALSE;
|
||||
}
|
||||
if (print == TRUE)
|
||||
{
|
||||
snprintf(token, sizeof(token), "Gas_Phase %d.\t%.350s",
|
||||
gas_phase_ptr->Get_n_user(), gas_phase_ptr->Get_description().c_str());
|
||||
dup_print(token, FALSE);
|
||||
}
|
||||
|
||||
/* Try to obtain a solution pointer */
|
||||
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, gas_phase_ptr->Get_n_solution()));
|
||||
prep();
|
||||
k_temp(use.Get_solution_ptr()->Get_tc(), use.Get_solution_ptr()->Get_patm());
|
||||
set(TRUE);
|
||||
converge = model();
|
||||
converge1 = check_residuals();
|
||||
if (converge == ERROR || converge1 == ERROR)
|
||||
{
|
||||
/* free_model_allocs(); */
|
||||
error_msg
|
||||
("Model failed to converge for initial gas phase calculation.",
|
||||
STOP);
|
||||
}
|
||||
use.Set_gas_phase_ptr(gas_phase_ptr);
|
||||
gas_phase_ptr->Set_total_p(0);
|
||||
gas_phase_ptr->Set_total_moles(0);
|
||||
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
|
||||
{
|
||||
cxxGasComp * gc_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
|
||||
int k;
|
||||
phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE);
|
||||
if (phase_ptr->in == TRUE)
|
||||
{
|
||||
lp = -phase_ptr->lk;
|
||||
for (rxn_ptr = &phase_ptr->rxn_x.token[0] + 1;
|
||||
rxn_ptr->s != NULL; rxn_ptr++)
|
||||
{
|
||||
lp += rxn_ptr->s->la * rxn_ptr->coef;
|
||||
}
|
||||
phase_ptr->p_soln_x = exp(lp * LOG_10);
|
||||
gas_phase_ptr->Set_total_p(gas_phase_ptr->Get_total_p() + phase_ptr->p_soln_x);
|
||||
phase_ptr->moles_x = phase_ptr->p_soln_x *
|
||||
gas_phase_ptr->Get_volume() / (R_LITER_ATM * tk_x);
|
||||
gc_ptr->Set_moles(phase_ptr->moles_x);
|
||||
gas_phase_ptr->Set_total_moles(gas_phase_ptr->Get_total_moles() + phase_ptr->moles_x);
|
||||
if (phase_ptr->p_c || phase_ptr->t_c)
|
||||
PR = true;
|
||||
}
|
||||
else
|
||||
{
|
||||
phase_ptr->moles_x = 0;
|
||||
}
|
||||
}
|
||||
if (fabs(gas_phase_ptr->Get_total_p() - use.Get_solution_ptr()->Get_patm()) > 5)
|
||||
{
|
||||
snprintf(token, sizeof(token),
|
||||
"WARNING: While initializing gas phase composition by equilibrating:\n%s (%.2f atm) %s (%.2f atm).\n%s.",
|
||||
" Gas phase pressure",
|
||||
(double) gas_phase_ptr->Get_total_p(),
|
||||
"is not equal to solution-pressure",
|
||||
(double) use.Get_solution_ptr()->Get_patm(),
|
||||
" Pressure effects on solubility may be incorrect");
|
||||
dup_print(token, FALSE);
|
||||
}
|
||||
|
||||
print_gas_phase();
|
||||
if (pr.user_print)
|
||||
print_user_print();
|
||||
if (PR /*&& use.Get_gas_phase_ptr()->total_p > 1.0*/)
|
||||
{
|
||||
std::ostringstream msg;
|
||||
msg << "\nWhile initializing gas phase composition by equilibrating:\n";
|
||||
msg << " Found definitions of gas critical temperature and pressure.\n";
|
||||
msg << " Going to use Peng-Robinson in subsequent calculations.\n";
|
||||
screen_msg(msg.str().c_str());
|
||||
output_msg(msg.str().c_str());
|
||||
log_msg(msg.str().c_str());
|
||||
}
|
||||
xgas_save(n_user);
|
||||
punch_all();
|
||||
/* free_model_allocs(); */
|
||||
}
|
||||
Utilities::Rxn_copies(Rxn_gas_phase_map, n_user, last);
|
||||
}
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
initial_surfaces(int print)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/*
|
||||
* Go through list of surfaces, make initial calculations
|
||||
* for any marked "new" that are defined to be in equilibrium with a
|
||||
* solution.
|
||||
*/
|
||||
int last, n_user, print1;
|
||||
|
||||
state = INITIAL_SURFACE;
|
||||
set_use();
|
||||
print1 = TRUE;
|
||||
|
||||
//std::map<int, cxxSurface>::iterator it = Rxn_surface_map.begin();
|
||||
//for ( ; it != Rxn_surface_map.end(); it++)
|
||||
//{
|
||||
//for (size_t nn = 0; nn < Rxn_new_surface.size(); nn++)
|
||||
//{
|
||||
//std::map<int, cxxSurface>::iterator it = Rxn_surface_map.find(Rxn_new_surface[nn]);
|
||||
for (std::set<int>::const_iterator nit = Rxn_new_surface.begin(); nit != Rxn_new_surface.end(); nit++)
|
||||
{
|
||||
std::map<int, cxxSurface>::iterator it = Rxn_surface_map.find(*nit);
|
||||
if (it == Rxn_surface_map.end())
|
||||
{
|
||||
assert(false);
|
||||
}
|
||||
cxxSurface * surface_ptr = &it->second;
|
||||
if (!surface_ptr->Get_new_def())
|
||||
continue;
|
||||
n_user = surface_ptr->Get_n_user();
|
||||
last = surface_ptr->Get_n_user_end();
|
||||
surface_ptr->Set_n_user_end(n_user);
|
||||
if (surface_ptr->Get_solution_equilibria())
|
||||
{
|
||||
if (print1 == TRUE && print == TRUE)
|
||||
{
|
||||
dup_print
|
||||
("Beginning of initial surface-composition calculations.",
|
||||
TRUE);
|
||||
print1 = FALSE;
|
||||
}
|
||||
if (print == TRUE)
|
||||
{
|
||||
std::ostringstream msg;
|
||||
msg << "Surface " << n_user << ".\t" << surface_ptr->Get_description().c_str();
|
||||
dup_print(msg.str().c_str(), FALSE);
|
||||
}
|
||||
use.Set_surface_ptr(surface_ptr);
|
||||
dl_type_x = use.Get_surface_ptr()->Get_dl_type();
|
||||
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, surface_ptr->Get_n_solution()));
|
||||
if (use.Get_solution_ptr() == NULL)
|
||||
{
|
||||
error_msg
|
||||
("Solution not found for initial surface calculation",
|
||||
STOP);
|
||||
}
|
||||
set_and_run_wrapper(-1, FALSE, FALSE, -1, 0.0);
|
||||
species_list_sort();
|
||||
print_surface();
|
||||
if (pr.user_print)
|
||||
print_user_print();
|
||||
punch_all();
|
||||
xsurface_save(n_user);
|
||||
/* free_model_allocs(); */
|
||||
}
|
||||
Utilities::Rxn_copies(Rxn_surface_map, n_user, last);
|
||||
}
|
||||
return (OK);
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::
|
||||
reactions(void)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
{
|
||||
/*
|
||||
* Make all reaction calculation which could include:
|
||||
* equilibrium with a pure-phase assemblage,
|
||||
* equilibrium with an exchanger,
|
||||
* equilibrium with an surface,
|
||||
* equilibrium with a gas phase,
|
||||
* equilibrium with a solid solution assemblage,
|
||||
* kinetics,
|
||||
* change of temperature,
|
||||
* mixture,
|
||||
* or irreversible reaction.
|
||||
*/
|
||||
int count_steps, use_mix;
|
||||
char token[2 * MAX_LENGTH];
|
||||
class save save_data;
|
||||
LDBLE kin_time;
|
||||
cxxKinetics *kinetics_ptr;
|
||||
|
||||
state = REACTION;
|
||||
/* last_model.force_prep = TRUE; */
|
||||
if (set_use() == FALSE)
|
||||
return (OK);
|
||||
/*
|
||||
* Find maximum number of steps
|
||||
*/
|
||||
dup_print("Beginning of batch-reaction calculations.", TRUE);
|
||||
count_steps = 1;
|
||||
if (use.Get_reaction_in() == TRUE && use.Get_reaction_ptr() != NULL)
|
||||
{
|
||||
cxxReaction *reaction_ptr = use.Get_reaction_ptr();
|
||||
if (reaction_ptr->Get_reaction_steps() > count_steps)
|
||||
count_steps = reaction_ptr->Get_reaction_steps();
|
||||
}
|
||||
if (use.Get_kinetics_in() == TRUE && use.Get_kinetics_ptr() != NULL)
|
||||
{
|
||||
if (use.Get_kinetics_ptr()->Get_reaction_steps() > count_steps)
|
||||
count_steps = use.Get_kinetics_ptr()->Get_reaction_steps();
|
||||
}
|
||||
if (use.Get_temperature_in() == TRUE && use.Get_temperature_ptr() != NULL)
|
||||
{
|
||||
int count = use.Get_temperature_ptr()->Get_countTemps();
|
||||
if (count > count_steps)
|
||||
{
|
||||
count_steps = count;
|
||||
}
|
||||
}
|
||||
if (use.Get_pressure_in() == TRUE && use.Get_pressure_ptr() != NULL)
|
||||
{
|
||||
int count = use.Get_pressure_ptr()->Get_count();
|
||||
if (count > count_steps)
|
||||
{
|
||||
count_steps = count;
|
||||
}
|
||||
}
|
||||
count_total_steps = count_steps;
|
||||
/*
|
||||
* save data for saving solutions
|
||||
*/
|
||||
// memcpy(&save_data, &save, sizeof(class save));
|
||||
save_data = save;
|
||||
/*
|
||||
*Copy everything to -2
|
||||
*/
|
||||
copy_use(-2);
|
||||
rate_sim_time_start = 0;
|
||||
rate_sim_time = 0;
|
||||
for (reaction_step = 1; reaction_step <= count_steps; reaction_step++)
|
||||
{
|
||||
overall_iterations = 0;
|
||||
snprintf(token, sizeof(token), "Reaction step %d.", reaction_step);
|
||||
if (reaction_step > 1 && incremental_reactions == FALSE)
|
||||
{
|
||||
copy_use(-2);
|
||||
}
|
||||
set_initial_moles(-2);
|
||||
dup_print(token, FALSE);
|
||||
/*
|
||||
* Determine time step for kinetics
|
||||
*/
|
||||
kin_time = 0.0;
|
||||
|
||||
if (use.Get_kinetics_in() == TRUE)
|
||||
{
|
||||
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, -2);
|
||||
kin_time = kinetics_ptr->Current_step((incremental_reactions==TRUE), reaction_step);
|
||||
}
|
||||
if (incremental_reactions == FALSE ||
|
||||
(incremental_reactions == TRUE && reaction_step == 1))
|
||||
{
|
||||
use_mix = TRUE;
|
||||
}
|
||||
else
|
||||
{
|
||||
use_mix = FALSE;
|
||||
}
|
||||
/*
|
||||
* Run reaction step
|
||||
*/
|
||||
run_reactions(-2, kin_time, use_mix, 1.0);
|
||||
if (incremental_reactions == TRUE)
|
||||
{
|
||||
rate_sim_time_start += kin_time;
|
||||
rate_sim_time = rate_sim_time_start;
|
||||
}
|
||||
else
|
||||
{
|
||||
rate_sim_time = kin_time;
|
||||
}
|
||||
if (state != ADVECTION)
|
||||
{
|
||||
punch_all();
|
||||
print_all();
|
||||
}
|
||||
/* saves back into -2 */
|
||||
if (reaction_step < count_steps)
|
||||
{
|
||||
saver();
|
||||
}
|
||||
}
|
||||
/*
|
||||
* save end of reaction
|
||||
*/
|
||||
// memcpy(&save, &save_data, sizeof(class save));
|
||||
save = save_data;
|
||||
if (use.Get_kinetics_in() == TRUE)
|
||||
{
|
||||
Utilities::Rxn_copy(Rxn_kinetics_map, -2, use.Get_n_kinetics_user());
|
||||
}
|
||||
saver();
|
||||
|
||||
/* free_model_allocs(); */
|
||||
//// set pr_in to false for following steps...
|
||||
// if (use.Get_pp_assemblage_in())
|
||||
// {
|
||||
// for (int i = 0; i < count_unknowns; i++)
|
||||
// {
|
||||
// if (x[i]->type == PP)
|
||||
// x[i]->phase->pr_in = false;
|
||||
// }
|
||||
// }
|
||||
// if (use.Get_gas_phase_in())
|
||||
// {
|
||||
// cxxGasPhase *gas_phase_ptr = use.Get_gas_phase_ptr();
|
||||
// for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
|
||||
// {
|
||||
// cxxGasComp *gc_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
|
||||
// int k;
|
||||
// class phase *phase_ptr = phase_bsearch(gc_ptr->Get_phase_name().c_str(), &k, FALSE);
|
||||
// assert(phase_ptr);
|
||||
// phase_ptr->pr_in = false;
|
||||
// }
|
||||
// }
|
||||
/* last_model.force_prep = TRUE; */
|
||||
rate_sim_time = 0;
|
||||
return (OK);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
int Phreeqc::initial_solutions(int print)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user