Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/reactmicp/micpsolvers/micpsolver.hpp
87604b01bf0dreactmicp/unsaturated

Last Change
Raw File

micpsolver.hpp
View Options

#ifndef SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
#define SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
#include <memory>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include "reactmicp/common.hpp"
#include "micpsolver_structs.hpp"
namespace specmicp {
namespace reactmicp {
namespace micpsolver {
//! MiCP solver for reactive transport
template <class program_t>
class MiCPParabolicSolver
{
public:
using program_ptr_t = std::shared_ptr<program_t>;
MiCPParabolicSolver(program_ptr_t the_program):
m_program(the_program)
{}
//! \brief Solve the program
MiCPParabolicSolverReturnCode solve_one_timestep(double dt,
ParabolicVariables& x);
int linesearch(Eigen::VectorXd& update, ParabolicVariables& x,
const Eigen::VectorXd &predictor, double dt);
MiCPParabolicSolverReturnCode compute_search_direction(Eigen::VectorXd& update,
ParabolicVariables &x,
double dt);
//! \brief Compute and reformulate the residuals
void setup_residuals(const ParabolicVariables& x);
//! \brief Compute the residuals
void compute_residuals(const ParabolicVariables& x,
Eigen::VectorXd& residual);
//! \brief Reformulate the residual
void reformulate_residuals(const ParabolicVariables& x);
//! \brief Reformulate the residual *in place*
void reformulate_residuals_inplace(const ParabolicVariables& x, Eigen::VectorXd& residual);
//! \brief Setup the jacobian
void setup_jacobian(ParabolicVariables& x,
const double dt,
Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian);
//! \brief Compute the jacobian
void compute_jacobian(ParabolicVariables& x,
const double dt,
Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian);
//! \brief Reformulate the jacobian
void reformulate_jacobian(const ParabolicVariables& x,
Eigen::SparseMatrix<double, Eigen::RowMajor> &jacobian, double dt);
MiCPParabolicSolverReturnCode check_convergence(int nb_iterations,
const Eigen::VectorXd& update,
const ParabolicVariables& solution,
const Eigen::VectorXd& residuals,
bool may_have_converged);
//! \brief Return the program
program_ptr_t get_program() {return m_program;}
//! \brief Return a reference to the options
MiCPParabolicSolverOptions& get_options() {return m_options;}
//! \brief Return a reference to the performance
MiCPParabolicSolverPerformance& get_performance() {return m_performance;}
//! \brief Scale the step if it is too long
double is_step_too_long(Eigen::VectorXd& update);
//! \brief Return the number of equations
double get_neq() {return m_program->get_neq();}
private:
//! \brief Compute the residuals in the linesearch
double compute_residuals_linesearch(const Eigen::VectorXd& update,
ParabolicVariables& vars,
const Eigen::VectorXd &predictor,
double dt,
double lambda);
program_ptr_t m_program;
MiCPParabolicSolverOptions m_options;
MiCPParabolicSolverPerformance m_performance;
Eigen::VectorXd m_residuals;
Eigen::VectorXd m_phi_residuals;
Eigen::VectorXd m_gradient;
Eigen::VectorXd m_predictor;
Eigen::VectorXd m_cp_variables;
std::vector<double> m_max_merits; //!< Contains the m best value of the merit function
bool m_gradient_step_taken;
double m_norm_0;
};
} // end namespace micpsolver
} // end namespace reactmicp
} // end namespace specmicp
// ============== //
// //
// implementation //
// //
// ============== //
#include "micpsolver.inl"
#endif // SPECMICP_REACTMICP_MICPSOLVER_MICPSOLVER_HPP
c4science · Help