Page MenuHomec4science
Files F81709619

smart_solver.cpp
No OneTemporary
Actions

File Metadata

Created
Sat, Sep 7, 13:28

smart_solver.cpp
View Options

/* =============================================================================
Copyright (c) 2014 - 2016
F. Georget <fabieng@princeton.edu> 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 "smart_solver.hpp"
#include "specmicp_database/data_container.hpp"
#include "specmicp/adimensional/adimensional_system_solution.hpp"
#include "specmicp/adimensional/adimensional_system_solver.hpp"
#include "specmicp/problem_solver/reactant_box.hpp"
#include "specmicp_common/log.hpp"
namespace specmicp {
struct SPECMICP_DLL_LOCAL SolverInitializer
{
bool has_been_set {false};
scalar_t vol_frac_water {-1};
std::unordered_map<std::string, scalar_t> log_molality_components;
std::unordered_map<std::string, scalar_t> vol_frac_minerals;
};
struct SmartAdimSolver::SmartAdimSolverImpl
{
std::shared_ptr<database::DataContainer> m_database;
AdimensionalSystemSolution m_solution;
AdimensionalSystemConstraints m_orig_constraints;
AdimensionalSystemSolverOptions m_options;
AdimensionalSystemConstraints m_current_constraints;
micpsolver::MiCPPerformance m_tot_perf;
SolverInitializer m_init_values;
bool m_is_solved {false};
SmartAdimSolverImpl(
std::shared_ptr<database::DataContainer> raw_data,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
):
m_database(raw_data),
m_orig_constraints(constraints),
m_options(options)
{
}
// Solver
// ------
// entry point of the algorithm
bool solve();
// all solve_* function will call this one
bool solve_with_constraints(const AdimensionalSystemConstraints& constraints);
// functions for all the different cases
bool solve_with_previous_solution();
bool solve_without_previous_solution();
bool solve_pressure_model();
bool solve_no_pressure_model();
bool solve_without_pressure_model();
bool solve_with_pressure_model();
bool solve_saturated();
// initializer
// ----------
//setter
void set_init_volfrac_water(scalar_t volume_fraction);
void set_init_molality(scalar_t molality);
void set_init_molality(std::unordered_map<std::string, scalar_t>&& molalities);
void set_init_volfrac_mineral(std::unordered_map<std::string, scalar_t>&& volume_fraction);
// initializer
void initialize_solution_vector(AdimensionalSystemSolver& solver, Vector& x);
};
SmartAdimSolver::SmartAdimSolver(
std::shared_ptr<database::DataContainer> raw_data,
const AdimensionalSystemConstraints& constraints,
const AdimensionalSystemSolverOptions& options
):
m_impl(utils::make_pimpl<SmartAdimSolverImpl>(raw_data, constraints, options))
{
}
SmartAdimSolver::SmartAdimSolver(
std::shared_ptr<database::DataContainer> raw_data,
const ReactantBox& reactant_box,
bool modify_db
):
SmartAdimSolver(
raw_data,
reactant_box.get_constraints(modify_db),
AdimensionalSystemSolverOptions()
)
{}
SmartAdimSolver::SmartAdimSolver(
std::shared_ptr<database::DataContainer> raw_data,
const ReactantBox& reactant_box,
const AdimensionalSystemSolverOptions& options,
bool modify_db
):
SmartAdimSolver(
raw_data,
reactant_box.get_constraints(modify_db),
options
)
{}
AdimensionalSystemSolverOptions& SmartAdimSolver::get_options() {
return m_impl->m_options;
}
SmartAdimSolver::~SmartAdimSolver() = default;
AdimensionalSystemSolution SmartAdimSolver::get_solution()
{
if (is_solved())
return m_impl->m_solution;
else
return AdimensionalSystemSolution();
}
bool SmartAdimSolver::is_solved()
{
return m_impl->m_is_solved;
}
bool SmartAdimSolver::solve()
{
return m_impl->solve();
}
void SmartAdimSolver::set_init_volfrac_water(scalar_t volume_fraction)
{
m_impl->set_init_volfrac_water(volume_fraction);
}
void SmartAdimSolver::set_init_molality(scalar_t molality)
{
m_impl->set_init_molality(molality);
}
void SmartAdimSolver::set_init_molality(std::unordered_map<std::string, scalar_t> molalities)
{
m_impl->set_init_molality(std::move(molalities));
}
void SmartAdimSolver::set_init_volfrac_mineral(std::unordered_map<std::string, scalar_t> volume_fraction)
{
m_impl->set_init_volfrac_mineral(std::move(volume_fraction));
}
void SmartAdimSolver::set_warmstart(const AdimensionalSystemSolution& previous_solution)
{
m_impl->m_solution = previous_solution;
}
void SmartAdimSolver::set_warmstart(AdimensionalSystemSolution&& previous_solution)
{
m_impl->m_solution = std::move(previous_solution);
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve()
{
DEBUG << "Entering smart solver";
if (m_solution.is_valid) {
if (solve_with_previous_solution())
{
return true;
}
else
{
ERROR << "Failed to solve with previous solution. "
<< "Try getting solution from scratch.";
}
}
return solve_without_previous_solution();
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_without_previous_solution()
{
if (m_orig_constraints.water_partial_pressure.use_partial_pressure_model)
{
return solve_pressure_model();
}
else
{
return solve_no_pressure_model();
}
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_with_previous_solution()
{
return solve_with_constraints(m_orig_constraints);
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_pressure_model()
{
bool loop_check = false; // variable to check if we entered a loop
restart:
if (solve_without_pressure_model())
{
WARNING << "Succeed to solve without pressure model";
m_options.allow_restart = false;
return solve_with_pressure_model();
}
else
{
ERROR << "Failed to solve problem without pressure model.";
if (loop_check)
{
m_is_solved = false;
return m_is_solved;
}
ERROR << "Smart solver, failed to solve problem without pressure model,"
" try to solve saturated problem first";
if (solve_saturated()) // try solving saturated as a seed
{
m_is_solved = false;
goto restart;
}
else // FAIL
{
ERROR << "Fail to solve system in saturated conditions";
m_is_solved = false;
return false;
}
}
return m_is_solved;
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_no_pressure_model()
{
// Try solving system as it is
if (solve_with_constraints(m_orig_constraints))
{
return true;
}
// if failed :
// -> if unsaturated, try solving saturated system first
// -> if saturated, give up
else if (m_orig_constraints.water_equation == WaterEquationType::MassConservation)
{
ERROR << "Error solving unsaturated system"
", try solving saturated system first.";
if (solve_saturated())
{
return solve_with_constraints(m_orig_constraints);
}
else
{
m_is_solved = false;
return false;
}
}
else
{
m_is_solved = false;
return false;
}
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_saturated()
{
if (m_orig_constraints.water_equation == WaterEquationType::MassConservation)
{
m_current_constraints = m_orig_constraints;
m_current_constraints.water_equation = WaterEquationType::SaturatedSystem;
m_current_constraints.water_partial_pressure.use_partial_pressure_model = false;
return solve_with_constraints(m_current_constraints);
}
return false;
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_without_pressure_model()
{
WARNING << "Temporarily remove pressure model";
m_current_constraints = m_orig_constraints;
m_current_constraints.water_partial_pressure.use_partial_pressure_model = false;
return solve_with_constraints(m_current_constraints);
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_with_pressure_model()
{
return solve_with_constraints(m_orig_constraints);
}
bool SmartAdimSolver::SmartAdimSolverImpl::solve_with_constraints(
const AdimensionalSystemConstraints& constraints
)
{
Vector x;
micpsolver::MiCPPerformance perf;
AdimensionalSystemSolver solver;
if (m_solution.is_valid)
{
solver = AdimensionalSystemSolver(m_database, constraints, m_solution, m_options);
x = m_solution.main_variables;
perf = solver.solve(x, false);
}
else
{
solver = AdimensionalSystemSolver(m_database, constraints, m_options);
if (m_init_values.has_been_set)
{
initialize_solution_vector(solver, x);
perf = solver.solve(x, false);
}
else
{
perf = solver.solve(x, true);
}
}
m_tot_perf += perf;
if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success)
{
m_is_solved = false;
}
else
{
m_solution = solver.get_raw_solution(x);
m_is_solved = true;
}
return m_is_solved;
}
void SmartAdimSolver::SmartAdimSolverImpl::set_init_volfrac_water(
scalar_t volume_fraction
)
{
m_init_values.has_been_set = true;
m_init_values.vol_frac_water = volume_fraction;
}
void SmartAdimSolver::SmartAdimSolverImpl::set_init_molality(scalar_t molality)
{
m_init_values.has_been_set = true;
scalar_t log_molality = std::log10(molality);
for (auto ind: m_database->range_aqueous_component())
{
m_init_values.log_molality_components[
m_database->get_label_component(ind)] = log_molality;
}
}
void SmartAdimSolver::SmartAdimSolverImpl::set_init_molality(
std::unordered_map<std::string, scalar_t>&& molalities
)
{
m_init_values.has_been_set = true;
for (auto elem: molalities)
{
m_init_values.log_molality_components[elem.first]
= std::log10(elem.second);
}
}
void SmartAdimSolver::SmartAdimSolverImpl::set_init_volfrac_mineral(
std::unordered_map<std::string, scalar_t>&& volume_fraction
)
{
m_init_values.has_been_set = true;
for (auto elem: volume_fraction)
{
m_init_values.vol_frac_minerals[elem.first] = elem.second;
}
}
void SmartAdimSolver::SmartAdimSolverImpl::initialize_solution_vector(
AdimensionalSystemSolver& solver,
Vector& x
)
{
// first : check the validity of the inputs
if (m_init_values.vol_frac_water < 0)
{
m_init_values.vol_frac_water
= m_options.system_options.restart_water_volume_fraction;
}
if (m_init_values.log_molality_components.empty())
{
set_init_molality(m_options.system_options.restart_concentration);
}
if (m_init_values.vol_frac_minerals.empty())
{
solver.initialize_variables(x,
m_init_values.vol_frac_water,
m_init_values.log_molality_components
);
}
else
{
solver.initialize_variables(x,
m_init_values.vol_frac_water,
m_init_values.log_molality_components,
m_init_values.vol_frac_minerals
);
}
}
} // end namespace specmicp

Event Timeline

c4science · Help