Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/dfpmsolver/parabolic_driver.hpp
522d1d188e48reactmicp/unsaturated

Last Change
Raw File

parabolic_driver.hpp
View Options

#ifndef SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#include "common.hpp"
#include "utils/sparse_solvers/sparse_solver.hpp"
#include "driver.hpp"
#include "parabolic_structs.hpp"
namespace specmicp {
namespace dfpmsolver {
//! \brief The parabolic driver
//!
//! Parabolic driver for finite element programs
template <class Program>
class ParabolicDriver: public Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>
{
public:
using base=Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>;
using base::program;
using base::get_neq;
using base::get_options;
using base::get_perfs;
using base::scaling;
using base::initialize_scaling;
ParabolicDriver(Program& the_program):
Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>(the_program),
m_velocity(Vector::Zero(the_program.get_tot_ndf())),
m_solver(nullptr)
{}
//! \brief Solve a timestep of length dt
ParabolicDriverReturnCode solve_timestep(scalar_t dt, Vector& displacement);
//! \brief Restart the current timestep
ParabolicDriverReturnCode restart_timestep(Vector& displacement);
//! \brief Check if the solution has been found
ParabolicDriverReturnCode check_convergence();
//! \brief Compute the residuals
void compute_residuals(const Vector& displacement,
Vector& residual)
{
compute_residuals(displacement, m_velocity, residual);
}
//! \brief Compute the residuals
void compute_residuals(const Vector& displacement,
const Vector& velocity,
Vector& residual)
{
program().compute_residuals(displacement, velocity, residual);
get_perfs().nb_call_residuals += 1;
}
//! \brief Compute the jacobian
void compute_jacobian(Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian
);
//! \brief Return the norm of the current residuals
double norm_residuals() {return m_residuals.norm();}
//! \brief Read/Write reference to the velocity vector
const Vector& get_velocity() const {return m_velocity;}
Vector& velocity() {return m_velocity;}
void set_velocity(Vector& velocity_vector);
//! \brief Reset the Sparse solver
//! This function should be called when a new solver need to be used
void reset_solver() {return m_solver.reset(nullptr);}
//! \brief Call this function if the pattern of the jacobian has changed
void jacobian_pattern_has_changed() {reset_solver();}
//! \brief Initialize the computation
void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement);
//! \brief Strang Linesearch
//!
//! ref :
//! - Matthies et al. (1979)
//! - JHP course notes
ParabolicLinesearchReturnCode strang_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda
);
scalar_t residuals_0() {return m_norm_0;}
private:
void reset_velocity();
//! \brief Backtracking Linesearch
//!
//! ref :
//! - Algo A6.3.1 : Dennis and Schnabel (1983)
//! - Nocedal & Wrigth (2006)
ParabolicLinesearchReturnCode backtracking_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda_out
);
//! Update the variables in displacement
void update_variable(
const Vector& update,
scalar_t lambda,
Vector& displacement
);
//! \brief Set the predictor
void set_predictor(Vector& displacement);
//! \brief perform the linesearch
ParabolicDriverReturnCode linesearch(
Vector &update,
Vector &displacements
);
//! Compute the residuals for the linesearch
double compute_residuals_linesearch(
Vector& update,
scalar_t lambda,
Vector& displacement
);
double compute_residuals_strang_linesearch(
Vector &update,
scalar_t lambda,
Vector &displacement
);
scalar_t update_norm(const Vector& update);
double m_norm_0;
double m_current_dt;
Eigen::VectorXd m_gradient;
Eigen::VectorXd m_velocity;
Eigen::VectorXd m_residuals;
Eigen::VectorXd m_predictor;
Eigen::SparseMatrix<double> m_jacobian;
sparse_solvers::SparseSolverPtr<Eigen::SparseMatrix<double>, Vector, Vector> m_solver;
bool m_velocity_is_initialized;
};
} // end namespace dfpmsolver
} // end namespace specmicp
// implementation
// ==============
#include "parabolic_driver.inl"
#endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
c4science · Help