From 643f53890e6fd18ed14b9c8602dc3581e127fdd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Max=20L=C3=BCbke?= Date: Tue, 28 Oct 2025 11:36:58 +0100 Subject: [PATCH] refactor(initialization): relocate initial setup methods --- src/phreeqcpp/mainsubs.cpp | 591 ------------------------------------- 1 file changed, 591 deletions(-) diff --git a/src/phreeqcpp/mainsubs.cpp b/src/phreeqcpp/mainsubs.cpp index 520c67db..09e3bdde 100644 --- a/src/phreeqcpp/mainsubs.cpp +++ b/src/phreeqcpp/mainsubs.cpp @@ -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::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::const_iterator nit = Rxn_new_solution.begin(); nit != Rxn_new_solution.end(); nit++) - { - std::map::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::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::iterator it = Rxn_exchange_map.find(Rxn_new_exchange[nn]); - for (std::set::const_iterator nit = Rxn_new_exchange.begin(); nit != Rxn_new_exchange.end(); nit++) - { - std::map::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::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::iterator it = Rxn_gas_phase_map.find(Rxn_new_gas_phase[nn]); - for (std::set::const_iterator nit = Rxn_new_gas_phase.begin(); nit != Rxn_new_gas_phase.end(); nit++) - { - std::map::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::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::iterator it = Rxn_surface_map.find(Rxn_new_surface[nn]); - for (std::set::const_iterator nit = Rxn_new_surface.begin(); nit != Rxn_new_surface.end(); nit++) - { - std::map::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) /* ---------------------------------------------------------------------- */