Page MenuHomec4science
Files F61994109

parabolic_driver.inl
No OneTemporary
Actions

File Metadata

Created
Fri, May 10, 06:45

parabolic_driver.inl
View Options

#include "parabolic_driver.hpp" // for syntaxic coloration only
#include "utils/log.hpp"
#include "sparse_solver.hpp"
namespace specmicp {
namespace dfpmsolver {
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::solve_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
initialize_timestep(dt, displacement);
compute_residuals(displacement, m_velocity, m_residuals);
m_norm_0 = norm_residuals();
get_perfs().nb_iterations = 0;
return restart_timestep(displacement);
}
template <class Program>
void ParabolicDriver<Program>::initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
m_residuals = Eigen::VectorXd::Zero(get_neq());
m_current_dt = dt;
set_predictor(displacement);
m_velocity.setZero();
program().apply_bc(dt, displacement, m_velocity);
}
template <class Program>
void ParabolicDriver<Program>::set_predictor(Vector& displacement)
{
if (get_options().alpha < 1)
m_predictor = displacement + (1-get_options().alpha)*m_current_dt*m_velocity;
else
m_predictor = displacement;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::restart_timestep(Vector& displacement)
{
ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet;
while (return_code == ParabolicDriverReturnCode::NotConvergedYet)
{
Eigen::VectorXd update(get_neq());
compute_residuals(displacement, m_velocity, m_residuals);
get_perfs().current_residual = m_residuals.norm();
return_code = check_convergence(-1.0);
if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break;
get_perfs().nb_iterations += 1;
compute_jacobian(displacement, m_velocity, m_jacobian);
m_gradient = m_jacobian.transpose()*m_residuals;
int retcode = solve_linear_system(update);
if (retcode <1)
{
ERROR << "Error when solving linear system : " << retcode << std::endl;
return_code = ParabolicDriverReturnCode::ErrorLinearSystem;
}
linesearch(update, displacement);
}
return return_code;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::check_convergence(scalar_t norm_update)
{
ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet;
const int nb_iterations = get_perfs().nb_iterations;
const scalar_t norm_residuals = m_residuals.norm();
DEBUG << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
if (norm_residuals/m_norm_0 < get_options().residuals_tolerance)
{
termcode = ParabolicDriverReturnCode::ResidualMinimized;
}
else if (nb_iterations > 0 and norm_update >0 and norm_update < get_options().step_tolerance)
{
if (norm_residuals > get_options().threshold_stationary_point)
{
ERROR << "Stationary point detected !";
termcode = ParabolicDriverReturnCode::StationaryPoint;
}
WARNING << "Error is minimized - may indicate a stationnary point";
termcode = ParabolicDriverReturnCode::ErrorMinimized;
}
else if (nb_iterations > get_options().maximum_iterations)
{
ERROR << "Maximum number of iteration reached (" << get_options().maximum_iterations << ")";
termcode = ParabolicDriverReturnCode::MaxIterations;
}
else if (get_perfs().maximum_step_taken)
{
get_perfs().nb_consecutive_max_step_taken += 1;
get_perfs().nb_max_step_taken += 1;
if (get_perfs().nb_consecutive_max_step_taken == get_options().max_iterations_at_max_length) {
ERROR << "Divergence detected - Maximum step length taken two many times";
termcode = ParabolicDriverReturnCode::MaxStepTakenTooManyTimes;
}
}
else
{
get_perfs().nb_consecutive_max_step_taken = 0;
}
get_perfs().return_code = termcode;
return termcode;
}
template <class Program>
int ParabolicDriver<Program>::solve_linear_system(Eigen::VectorXd& update)
{
return solve_sparse_linear_system(m_residuals, m_jacobian, update, get_options().sparse_solver);
}
template <class Program>
double ParabolicDriver<Program>::compute_residuals_linesearch(
Vector &update,
scalar_t lambda,
Vector &displacement
)
{
Eigen::VectorXd velocity(m_velocity);
Eigen::VectorXd residual = Eigen::VectorXd::Zero(get_neq());
program().update_solution(update, lambda, get_options().alpha*m_current_dt, m_predictor, displacement, velocity);
compute_residuals(displacement, velocity, residual);
return 0.5*residual.squaredNorm();
}
template <class Program>
void ParabolicDriver<Program>::update_variable(
const Vector& update,
scalar_t lambda,
Vector& displacement
)
{
program().update_solution(update, lambda, get_options().alpha*m_current_dt,
m_predictor, displacement, m_velocity);
}
template <class Program>
int ParabolicDriver<Program>::linesearch(
Vector& update,
Vector& displacements
)
{
// References
// ----------
// - Algo A6.3.1 : Dennis and Schnabel (1983)
// - Nocedal & Wrigth (2006)
DEBUG << "Linesearch";
Eigen::VectorXd xp(displacements);
double fcp;
get_perfs().maximum_step_taken = false;
int retcode = 2;
const scalar_t alpha = 1e-4;
scalar_t newtlen = base::is_step_too_long(update);
scalar_t init_slope = m_gradient.dot(update);
const scalar_t rellength = update.cwiseAbs().maxCoeff();
const scalar_t minlambda = get_options().step_tolerance / rellength;
scalar_t lambda = 1.0;
scalar_t lambda_prev = lambda;
scalar_t merit_value = 0.5*m_residuals.squaredNorm();
// new residual
fcp = compute_residuals_linesearch(update, lambda, xp);
// Skip linesearch if enough progress is done
// ------------------------------------------
if (fcp < get_options().coeff_accept_newton_step *merit_value)
{
DEBUG << "Skip linesearch ";
update_variable(update, lambda, displacements);
get_perfs().current_update = lambda * update.norm();
return 0;
}
// The linesearch
// --------------
scalar_t fc = merit_value;
scalar_t fcp_prev;
int cnt = 0;
do
{
SPAM << "cnt : " <<cnt << " - lambda : " << lambda << " # " << fc << " # " << fcp << " # " << init_slope;
if (fcp <= fc + alpha*lambda*init_slope) //pg760 Fachinei2003
{
retcode = 0;
if (lambda ==1 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
break;
}
else if (lambda < minlambda)
{
retcode = 1;
break;
}
else
{
//WARNING << "denom : " << fcp - fc -init_slope;
// Select a new step length
// - - - - - - - - - - - -
double lambdatmp;
if (cnt == 0) { // only a quadratic at the first
lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
}
else
{
const scalar_t factor = 1.0 /(lambda - lambda_prev);
const scalar_t x1 = fcp - fc - lambda*init_slope;
const scalar_t x2 = fcp_prev - fc - lambda_prev*init_slope;
const scalar_t a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
const scalar_t b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
if (a == 0)
{ // cubic interpolation is in fact a quadratic interpolation
lambdatmp = - init_slope/(2*b);
}
else
{
const scalar_t disc = b*b-3*a*init_slope;
lambdatmp = (-b+std::sqrt(disc))/(3*a);
}
if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
}
//WARNING << "lambdatmp : " << lambdatmp;
lambda_prev = lambda;
fcp_prev = fcp;
if (lambdatmp < 0.1*lambda) {
lambda = 0.1 * lambda;
} else {
lambda = lambdatmp;
}
}
if (not std::isfinite(lambda))
{
ERROR << "Lambda is non finite - we stop";
throw std::runtime_error("Lambda is nonfinite : "+std::to_string(lambda));
}
fcp = compute_residuals_linesearch(update, lambda, xp);
//xp.velocity = x.velocity;
++cnt;
} while(retcode == 2 and cnt < 100);
//WARNING << "Lambda : " << lambda << " - iterations : " << cnt;
if (cnt == 100)
{
ERROR << "Too much linesearch iterations ! We stop";
}
get_perfs().current_update = lambda * update.norm();
update_variable(update, lambda, displacements);
return retcode;
}
} // end namespace dfpmsolver
} // end namespace specmicp

Event Timeline

c4science · Help