Page MenuHomec4science

parabolic_driver.hpp
No OneTemporary

File Metadata

Created
Tue, May 21, 18:43

parabolic_driver.hpp

#ifndef SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#define SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP
#include "common.hpp"
#include <Eigen/Sparse>
#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;
ParabolicDriver(Program& the_program):
Driver<Program, ParabolicDriverOptions, ParabolicDriverPerformance>(the_program)
{}
//! \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(scalar_t norm_update);
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
)
{
program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt);
get_perfs().nb_call_jacobian += 1;
}
//! \brief Return the norm of the current residuals
double norm_residuals() {return m_residuals.norm();}
Eigen::VectorXd& get_velocity() {return m_velocity;}
private:
void update_variable(
const Vector& update,
scalar_t lambda,
Vector& displacement
);
//! \brief Set the predictor
void set_predictor(Vector& displacement);
//! \brief solve the linear system
int solve_linear_system(Eigen::VectorXd& update);
//! \brief perform the linesearch
int linesearch(Vector &update,
Vector &displacements);
double compute_residuals_linesearch(
Vector& update,
scalar_t lambda,
Vector& displacement
);
//! \brief Initialize the computation
void initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement);
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;
};
} // end namespace dfpmsolver
} // end namespace specmicp
// implementation
// ==============
#include "parabolic_driver.inl"
#endif // SPECMICP_DFPMSOLVER_PARABOLICDRIVER_HPP

Event Timeline