diff --git a/src/specmicp/adimensional/adimensional_system_solver.cpp b/src/specmicp/adimensional/adimensional_system_solver.cpp index 80835b9..cc3d463 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.cpp +++ b/src/specmicp/adimensional/adimensional_system_solver.cpp @@ -1,522 +1,547 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #include "adimensional_system_solver.hpp" #include "adimensional_system.hpp" #include "adimensional_system_solution.hpp" #include "adimensional_system_pcfm.hpp" #include "specmicp_common/micpsolver/micpsolver.hpp" #include "specmicp_common/log.hpp" #include namespace specmicp { static constexpr int max_solver_restart = 3; // max number of time the solver is // restarted in case of failure // solve_equilibrium function // ========================== AdimensionalSystemSolution solve_equilibrium( std::shared_ptr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data, constraints, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } AdimensionalSystemSolution solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ) { AdimensionalSystemSolver solver(data_ptr, constraints, previous_solution, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) throw std::runtime_error("Failed to solve the problem"); return solver.get_raw_solution(variables); } // The solver // ========== // constructor // ----------- AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ): OptionsHandler(), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, get_options().system_options, get_options().units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolver::AdimensionalSystemSolver( RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ): OptionsHandler(options), m_data(data), m_system(std::make_shared( data, constraints, previous_solution, options.system_options, options.units_set )), m_var(Vector::Zero(data->nb_component()+data->nb_ssites()+1+data->nb_mineral())) {} AdimensionalSystemSolution AdimensionalSystemSolver::get_raw_solution( Vector& x ) { set_true_variable_vector(x); return m_system->get_solution(x, m_var); } AdimensionalSystemSolution AdimensionalSystemSolver::unsafe_get_raw_solution( Vector& x ) { set_true_variable_vector(x); return m_system->unsafe_get_solution(x, m_var); } +AdimensionalSystemSolution AdimensionalSystemSolver::merge_raw_solution_with_immobile_species( + Vector& x, + const AdimensionalSystemSolution& sol_immobile + ) +{ + set_true_variable_vector(x); + AdimensionalSystemSolution sol = m_system->unsafe_get_solution(x, m_var); + + sol.sorbed_molalities = sol_immobile.sorbed_molalities; + if (m_data->nb_mineral() > 0) + { + auto m_start = m_system->dof_mineral(0); + sol.main_variables.segment(m_start, m_data->nb_mineral()) = + sol_immobile.main_variables.segment(m_start, m_data->nb_mineral()); + } + if (m_data->nb_ssites() > 0) + { + auto s_start = m_system->dof_surface(0); + sol.main_variables.segment(s_start, m_data->nb_ssites()) = + sol_immobile.main_variables.segment(s_start, m_data->nb_ssites()); + + } + return sol; +} + // Solving the system // ------------------ micpsolver::MiCPPerformance AdimensionalSystemSolver::solve(Vector& x, bool init) { // reset the units // - - - - - - - - m_system->set_units(get_options().units_set); // initialize the system // - - - - - - - - - - - if (init) { m_system->reasonable_starting_guess(x); if (get_options().use_pcfm) { run_pcfm(x); } } else if (get_options().force_pcfm) { run_pcfm(x); } // Solve the system // - - - - - - - - micpsolver::MiCPPerformance perf = solve_system(x); // If failed : try again // - - - - - - - - - - - int cnt = 0; while (perf.return_code < micpsolver::MiCPSolverReturnCode::Success and get_options().allow_restart and cnt < max_solver_restart ) { WARNING << "Failed to solve the system ! Return code :" << (int) perf.return_code << ". We shake it up and start again"; // Try safer options const scalar_t save_penalization_factor = get_options().solver_options.penalization_factor; const bool save_scaling = get_options().solver_options.use_scaling; get_options().solver_options.use_scaling = true; if (save_penalization_factor == 1.0) get_options().solver_options.penalization_factor = 0.8; // Re-initialize the system set_return_vector(x); m_system->reasonable_restarting_guess(x); if (get_options().use_pcfm or get_options().force_pcfm) run_pcfm(x); // Solve the system micpsolver::MiCPPerformance perf2 = solve_system(x); // Restore options get_options().solver_options.penalization_factor = save_penalization_factor; get_options().solver_options.use_scaling = save_scaling; // Record the work performed perf += perf2; ++cnt; } // If solution is good // - - - - - - - - - - if (perf.return_code >= micpsolver::MiCPSolverReturnCode::Success) { set_return_vector(x); } return perf; } micpsolver::MiCPPerformance AdimensionalSystemSolver::solve_system(Vector& x) { set_true_variable_vector(x); if (not m_system->system_is_box_constrained()) { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_perfs(); } else { micpsolver::MiCPSolver solver(m_system); solver.set_options(get_options().solver_options); solver.solve(m_var); return solver.get_perfs(); } } // Variables management // --------------------- void AdimensionalSystemSolver::set_true_variable_vector(const Vector& x) { // This function sets the true variable for the specmicp system // the order is important // this is synchronised with the specmicp system, but needs to be set here const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i {0}; if (m_system->is_active_component(0)) { m_var(new_i) = x(m_system->dof_water()); ++new_i; } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(), i ); if (it != non_active_component.cend()) continue; scalar_t value = x(m_system->dof_component(i)); if (value == -HUGE_VAL) // check for previously undefined value { value = m_system->get_options().new_component_concentration; } m_var(new_i) = value; ++new_i; } if (m_system->surface_model() != SurfaceEquationType::NoEquation) { for (index_t q: m_data->range_ssites()) { if (m_system->ideq_surf(q) != no_equation) { m_var(new_i) = x(m_system->dof_surface(q)); ++new_i; } } if (m_system->ideq_surf_pot() != no_equation) { m_var(new_i) = x(m_system->dof_surface_potential()); ++new_i; } } if (m_system->ideq_electron() != no_equation) { m_var(new_i) = x(m_system->dof_electron()); ++new_i; } for (index_t m: m_data->range_mineral()) { if (m_system->is_active_mineral(m)) { m_var(new_i) = x(m_system->dof_mineral(m)); ++new_i; } } m_var.conservativeResize(new_i); specmicp_assert(new_i == (unsigned) m_system->total_variables()); } void AdimensionalSystemSolver::set_return_vector(Vector& x) { // This function sets the variable from the specmicp "raw" solution // the order is important // this is synchronised with the specmicp system, but needs to be set here const std::vector& non_active_component = m_system->get_non_active_component(); uindex_t new_i = 0; if (m_system->is_active_component(0)) { x(m_system->dof_water()) = m_var(new_i); ++new_i; } else { x(m_system->dof_water()) = m_system->volume_fraction_water(x); } for (index_t i: m_data->range_aqueous_component()) { const auto it = std::find(non_active_component.cbegin(), non_active_component.cend(),i); if (it != non_active_component.cend()) { x(m_system->dof_component(i)) = -HUGE_VAL; continue; } x(m_system->dof_component(i)) = m_var(new_i) ; ++new_i; } for (index_t q: m_data->range_ssites()) { if (m_system->ideq_surf(q) != no_equation) { x(m_system->dof_surface(q)) = m_var(new_i); ++new_i; } else { x(m_system->dof_surface(q)) = -HUGE_VAL; } } if (m_system->ideq_surf_pot() != no_equation) { x(m_system->dof_surface_potential()) = m_var(new_i); ++new_i; } else { x(m_system->dof_surface_potential()) = -HUGE_VAL; } if (m_system->ideq_electron() != no_equation) { x(m_system->dof_electron()) = m_var(new_i); ++new_i; } else x(m_system->dof_electron()) = -HUGE_VAL; for (index_t m: m_data->range_mineral()) { if (m_system->is_active_mineral(m)) { x(m_system->dof_mineral(m)) =m_var(new_i); ++new_i; } else { x(m_system->dof_mineral(m)) = 0.0; } } } // PCFM // ---- void AdimensionalSystemSolver::run_pcfm(Vector &x) { DEBUG << "Start PCFM initialization."; // we set up the true variable set_true_variable_vector(x); // The residual is computed to have some point of comparison Vector residuals(m_system->total_variables()); residuals.setZero(); m_system->get_residuals(m_var, residuals); const scalar_t res_0 = residuals.norm(); // the pcfm iterations are executed AdimensionalSystemPCFM pcfm_solver(m_system); PCFMReturnCode retcode = pcfm_solver.solve(m_var); // Check the answer if (retcode < PCFMReturnCode::Success) { // small prograss is still good enough m_system->get_residuals(m_var, residuals); const scalar_t final_residual = residuals.norm(); DEBUG << "Final pcfm residuals : " << final_residual << " set_secondary_variables(m_var); } } // Initialisation of variables // --------------------------- void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, std::unordered_map log_molalities, std::unordered_map volume_fraction_minerals, scalar_t log_free_sorption_site_concentration ) { m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; for (auto pair: log_molalities) { index_t idc = m_data->get_id_component(pair.first); if (idc == no_species or idc == m_data->electron_index() or idc == m_data->water_index()) { throw std::invalid_argument("This is not an aqueous component : "+pair.first); } if (pair.second > 0) { WARNING << "Initial molality for : " << pair.first << "is bigger than 1 molal."; } x(m_system->dof_component(idc)) = pair.second; } for (auto pair: volume_fraction_minerals) { index_t idm = m_data->get_id_mineral(pair.first); if (idm == no_species ) { throw std::invalid_argument("This is not a mineral at equilibrium : "+pair.first); } if (pair.second < 0 or pair.second > 1) { WARNING << "Initial volume fraction for : " << pair.first << "is not between 0 and 1."; } x(m_system->dof_mineral(idm)) = pair.second; } if (log_free_sorption_site_concentration != 0.0) { for (index_t q: m_data->range_ssites()) { x(m_system->dof_surface(q)) = log_free_sorption_site_concentration; } } } void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities ) { m_system->reasonable_starting_guess(x); if (volume_fraction_water < 0 or volume_fraction_water > 1) { WARNING << "Initial guess for the volume fraction of water is not between 0 and 1"; } x(m_system->dof_water()) = volume_fraction_water; if (log_molalities > 0) { WARNING << "Initial molality for : " << log_molalities << "is bigger than 1 molal."; } x.segment(1, m_data->nb_component()-1).setConstant(log_molalities); } void AdimensionalSystemSolver::initialize_variables( Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities, scalar_t log_free_sorption_site_concentration ) { initialize_variables(x, volume_fraction_water, log_molalities); if (log_free_sorption_site_concentration != 0.0) { for (index_t q: m_data->range_ssites()) { x(m_system->dof_surface(q)) = log_free_sorption_site_concentration; } } } } // end namespace specmicp diff --git a/src/specmicp/adimensional/adimensional_system_solver.hpp b/src/specmicp/adimensional/adimensional_system_solver.hpp index 1b4cfde..636b0e0 100644 --- a/src/specmicp/adimensional/adimensional_system_solver.hpp +++ b/src/specmicp/adimensional/adimensional_system_solver.hpp @@ -1,230 +1,238 @@ /* ============================================================================= Copyright (c) 2014 - 2016 F. Georget Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * ============================================================================= */ #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP #include "adimensional_system_solver_structs.hpp" #include "specmicp_database/database.hpp" #include "specmicp_common/options_handler.hpp" #include #include /*! \file adimensional_system_solver.hpp \brief Solve a reduced system \code{cpp} std::shared_ptr raw_data; AdimensionalSystemConstraints constraints; // ... AdimensionalSystemSolution previous_solution; // ... AdimensionalSystemSolverOptions options; // ... AdimensionalSystemSolver solver(raw_data, constraints, previous_solution, options); Vector variables; micpsolver::MiCPPerformance perf = solver.solve(variables, true); if (perf.return_code <= micpsolver::MiCPSolverReturnCode::NotConvergedYet) { throw std::runtime_error("Failed to solve the problem"); } AdimensionalSystemSolution solution = solver.get_raw_solution(variables); \endcode */ namespace specmicp { // forward declaration class AdimensionalSystem; struct AdimensionalSystemSolution; //! \brief Solve an adimensional system //! //! Take care of non-existing component in the system //! and restart the computation if necessary //! //! \ingroup specmicp_api class SPECMICP_DLL_PUBLIC AdimensionalSystemSolver: public OptionsHandler { public: //! \brief Pointer to the system using SystemPtr = std::shared_ptr; //! Default constructor //! //! Another constructor is necessary AdimensionalSystemSolver() {} //! \brief Initialise a solver from scratch //! //! \param data A raw database //! \param constraints The system to solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints ); //! \brief Initialise a solver from scratch //! //! \param data A raw database //! \param constraints The system to solver //! \param options (optional) customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database //! \param constraints The system to solver //! \param previous_solution A solution of a similar system that will be used to initialize the system AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution ); //! \brief Initialise a solver from a previous solution //! //! \param data A raw database //! \param constraints The system to solver //! \param previous_solution A solution of a similar system that will be used to initialize the system //! \param options customize the behavior of the solver AdimensionalSystemSolver(RawDatabasePtr data, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); //! \brief solve the problem using initial guess x //! //! \param[in,out] x in -> initial guess, out -> solution //! \param init if true, the algorithm guess a starting point micpsolver::MiCPPerformance solve(Vector& x, bool init=false); //! \brief Return the system used for the computation SystemPtr get_system() {return m_system;} //! \brief Return the solution in a manageable form //! //! \param x The solution (complete set of the main variables) AdimensionalSystemSolution get_raw_solution(Vector& x); //! \brief Return the solution in a manageable form //! //! \param x The solution (complete set of the main variables) AdimensionalSystemSolution unsafe_get_raw_solution(Vector& x); + //! \brief Return the solution merge with solid phases of another + //! + //! This is useful when used with the "immobilized solid phases" + AdimensionalSystemSolution merge_raw_solution_with_immobile_species( + Vector& x, + const AdimensionalSystemSolution& sol_immobile + ); + //! \brief Initialize the problem using the Positive continuous fraction method //! //! \sa specmicp::AdimensionalSystemPCFM void run_pcfm(Vector& x); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for chosen aqueous component //! \param volume_fraction_minerals volume fraction of the minerals //! \param log_free_sorption_site_concentration concentration of free sorption site void initialize_variables(Vector& x, scalar_t volume_fraction_water, std::unordered_map log_molalities, std::unordered_map volume_fraction_minerals = std::unordered_map(), scalar_t log_free_sorption_site_concentration = 0 ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities ); //! \brief Custom initialisation of variables //! //! The amount compon //! \param[out] x The initial guess (complete set of the main variables) //! \param volume_fraction_water volume fraction of water //! \param log_molalities log_10 of the molalities for all aqueous components //! \param log_free_sorption_site_concentration concentration of free sorption site void initialize_variables(Vector& x, scalar_t volume_fraction_water, scalar_t log_molalities, scalar_t log_free_sorption_site_concentration ); private: //! \brief set up the true variable vector //! //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_true_variable_vector(const Vector& x); //! \brief set up the true solution vector //! //! add zero components //! \param x The solution (complete set of the main variables) void SPECMICP_DLL_LOCAL set_return_vector(Vector& x); //! \brief solve the problem micpsolver::MiCPPerformance SPECMICP_DLL_LOCAL solve_system(Vector& x); RawDatabasePtr m_data; //! The raw database SystemPtr m_system; //! The system to solve Vector m_var; //! Copy of the solution vector (necessary in case of failing) }; //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints Constraints applied to the system //! \param options Options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolverOptions& options ); //! \brief Solve a reduced system, function provided for convenience //! //! \param data_ptr the database //! \param constraints constraints applied to the system //! \param previous_solution a previous solution //! \param options options for the solver AdimensionalSystemSolution SPECMICP_DLL_PUBLIC solve_equilibrium( RawDatabasePtr data_ptr, const AdimensionalSystemConstraints& constraints, const AdimensionalSystemSolution& previous_solution, const AdimensionalSystemSolverOptions& options ); } // end namespace specmicp #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLVER_HPP