Homec4science
Diffusion SpecMiCP / ReactMiCP (oxydoreduc) Browse

specmicp/src/dfpmsolver/parabolic_driver.inl
d4f8837fe33aoxydoreduc

Last Change
Raw File

parabolic_driver.inl
View Options

#include "parabolic_driver.hpp" // for syntaxic coloration only
#include "utils/log.hpp"
//#include "utils/sparse_solver.hpp"
namespace specmicp {
namespace dfpmsolver {
template <class Program>
void ParabolicDriver<Program>::compute_jacobian(Vector& displacement,
Vector& velocity,
Eigen::SparseMatrix<scalar_t>& jacobian
)
{
program().compute_jacobian(displacement, velocity, jacobian, get_options().alpha*m_current_dt);
jacobian = jacobian*scaling().asDiagonal();
jacobian.makeCompressed();
get_perfs().nb_call_jacobian += 1;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::solve_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
initialize_timestep(dt, displacement);
return restart_timestep(displacement);
}
template <class Program>
void ParabolicDriver<Program>::initialize_timestep(scalar_t dt, Eigen::VectorXd& displacement)
{
initialize_scaling();
m_residuals = Eigen::VectorXd::Zero(get_neq());
m_current_dt = dt;
set_predictor(displacement);
reset_velocity();
program().apply_bc(dt, displacement, m_velocity);
compute_residuals(displacement, m_velocity, m_residuals);
m_norm_0 = norm_residuals();
get_perfs().nb_iterations = 0;
}
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>
scalar_t ParabolicDriver<Program>::update_norm(const Vector& update)
{
// l-∞ scaled norm
scalar_t norm = 0.0;
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
const index_t id_eq = program().id_equation(dof);
if (id_eq == no_species or update(id_eq) == 0.0) continue;
norm = std::max(norm, std::abs(update(id_eq))/(std::max(std::abs(m_velocity(dof)), 1.0/scaling(id_eq))));
}
return norm;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::restart_timestep(Vector& displacement)
{
ParabolicDriverReturnCode return_code = ParabolicDriverReturnCode::NotConvergedYet;
//m_solver.reset(nullptr);
Eigen::VectorXd update(get_neq());
update.setZero();
get_perfs().current_update = 0;
bool force_recompute_jacobian = true;
while (return_code == ParabolicDriverReturnCode::NotConvergedYet)
{
compute_residuals(displacement, m_velocity, m_residuals);
get_perfs().absolute_residual = m_residuals.norm();
get_perfs().current_residual = m_residuals.norm()/m_norm_0;
get_perfs().current_update = update_norm(update);
DEBUG << " NB iterations : " << get_perfs().nb_iterations
<< " - res : " << get_perfs().current_residual
<< " - update : " << get_perfs().current_update;
return_code = check_convergence();
if (return_code != ParabolicDriverReturnCode::NotConvergedYet) break;
if (m_solver == nullptr)
{
m_solver = sparse_solvers::get_sparse_solver<
Eigen::SparseMatrix<scalar_t>, Vector, Vector>(get_options().sparse_solver);
compute_jacobian(displacement, m_velocity, m_jacobian);
m_solver->analyse_pattern(m_jacobian);
m_solver->decompose(m_jacobian);
force_recompute_jacobian =false;
}
else if (force_recompute_jacobian
or get_perfs().nb_iterations % get_options().quasi_newton == 0)
{
compute_jacobian(displacement, m_velocity, m_jacobian);
m_solver->decompose(m_jacobian);
}
get_perfs().nb_iterations += 1;
m_gradient = m_jacobian.transpose()*m_residuals;
sparse_solvers::SparseSolverReturnCode retcode = m_solver->solve_scaling(m_residuals, scaling(), update);
if (retcode != sparse_solvers::SparseSolverReturnCode::Success)
{
ERROR << "Error when solving linear system : " << (int) retcode << std::endl;
return ParabolicDriverReturnCode::ErrorLinearSystem;
}
//if (update.norm() < get_options().step_tolerance) return_code = ParabolicDriverReturnCode::ErrorMinimized;
//else
return_code = linesearch(update, displacement);
}
return return_code;
}
template <class Program>
ParabolicDriverReturnCode ParabolicDriver<Program>::check_convergence()
{
ParabolicDriverReturnCode termcode = ParabolicDriverReturnCode::NotConvergedYet;
const int nb_iterations = get_perfs().nb_iterations;
const scalar_t norm_residuals = get_perfs().current_residual;
const scalar_t norm_update = get_perfs().current_update;
//std::cout << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0 << std::endl;
DEBUG << "Residuals : " << nb_iterations << " - " << norm_residuals/m_norm_0;
if (norm_residuals < get_options().residuals_tolerance)
{
termcode = ParabolicDriverReturnCode::ResidualMinimized;
}
else if (get_perfs().absolute_residual < get_options().absolute_tolerance)
{
termcode = ParabolicDriverReturnCode::ResidualMinimized;
}
else if (nb_iterations > 0 and norm_update > 0.0 and norm_update < 1.01*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>
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>
double ParabolicDriver<Program>::compute_residuals_strang_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 update.dot(residual);
}
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>
ParabolicDriverReturnCode ParabolicDriver<Program>::linesearch(
Vector &update,
Vector &displacements
)
{
base::is_step_too_long(update);
get_perfs().maximum_step_taken = false;
scalar_t lambda;
ParabolicLinesearchReturnCode retcode;
switch (get_options().linesearch) {
case ParabolicLinesearch::Bactracking:
retcode = backtracking_linesearch(update, displacements, lambda);
break;
case ParabolicLinesearch::Strang:
retcode = strang_linesearch(update, displacements, lambda);
break;
default:
throw std::runtime_error("Linesearch type for Parabolic driver is not recognized");
break;
}
if (retcode != ParabolicLinesearchReturnCode::Success)
{
return ParabolicDriverReturnCode::LinesearchFailed;
}
update_variable(update, lambda, displacements);
update *= lambda;
return ParabolicDriverReturnCode::NotConvergedYet;
}
namespace internal
{
//! \brief Return true if a and b have the same sign
inline bool have_same_sign(double a, double b)
{
return (std::copysign(1., a) * std::copysign(1., b) > 0);
}
} // end namespace internal
template <class Program>
ParabolicLinesearchReturnCode ParabolicDriver<Program>::strang_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda
)
{
DEBUG << "Strang linesearch";
Eigen::VectorXd xp(displacements);
const scalar_t s_tol = 0.5;
const scalar_t s_max = 2.0;
const int lin_max = 10;
scalar_t s_b = 0.0;
scalar_t g_b = 0.5*m_residuals.squaredNorm();
scalar_t s_a = 1.0;
scalar_t g_a = compute_residuals_strang_linesearch(update, 1.0, xp);
scalar_t newtlen = (scaling().asDiagonal()*update).norm();
// const scalar_t s_r = s_a;
const scalar_t g_r = g_a;
if (std::abs(g_a) <= s_tol*std::abs(g_b))
{
DEBUG << "Skip linesearch ";
lambda = 1.0;
if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
return ParabolicLinesearchReturnCode::Success;
}
while (internal::have_same_sign(g_a, g_b) and s_a < s_max)
{
s_b = s_a;
s_a = 2*s_a;
g_b = g_a;
g_a = compute_residuals_strang_linesearch(update, s_a, xp);
}
scalar_t g_ = g_a;
scalar_t g_0 = g_a;
scalar_t s = s_a;
int l;
for (l=0; l<lin_max; ++l)
{
if (internal::have_same_sign(g_a, g_b) or
std::abs(g_) < s_tol*std::abs(g_0)
or std::abs(s_a - s_b) < s_tol*(s_a + s_b)/2 )
break;
s = s_a - g_a * ( s_a - s_b)/(g_a - g_b);
g_ = compute_residuals_strang_linesearch(update, s, xp);
if (internal::have_same_sign(g_, g_a))
g_b = g_b/2.;
else
{
s_b = s_a;
g_b = g_a;
}
s_a = s;
g_a = g_;
}
if (l == lin_max) return ParabolicLinesearchReturnCode::MaximumIterations;
if (g_a >= g_r) {
WARNING << "Failed to find better update in Strang linesearch";
//lambda = 0.1*s_r;
return backtracking_linesearch(update, displacements, lambda);
}
lambda = s_a;
if (lambda == 1.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
return ParabolicLinesearchReturnCode::Success;
}
template <class Program>
ParabolicLinesearchReturnCode ParabolicDriver<Program>::backtracking_linesearch(
Vector& update,
Vector& displacements,
scalar_t& lambda_out
)
{
// References
// ----------
// - Algo A6.3.1 : Dennis and Schnabel (1983)
// - Nocedal & Wrigth (2006)
DEBUG << "Linesearch";
Eigen::VectorXd xp(displacements);
double fcp;
int retcode = 2; // 2 not converged, 1 problem, 0 success
const scalar_t alpha = 1e-4;
scalar_t newtlen = (scaling().asDiagonal()*update).norm();
scalar_t init_slope = m_gradient.dot(update);
const scalar_t rellength = update_norm(update);
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 ";
lambda_out = 1.0;
return ParabolicLinesearchReturnCode::Success;
}
// 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.0 and (newtlen > 0.99 * get_options().maximum_step_length)) {
get_perfs().maximum_step_taken = true;
}
break;
}
else if (lambda < minlambda)
{
retcode = 0;
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";
return ParabolicLinesearchReturnCode::Divergence;
}
fcp = compute_residuals_linesearch(update, lambda, xp);
//xp.velocity = x.velocity;
++cnt;
} while(retcode == 2 and cnt < 50);
//WARNING << "Lambda : " << lambda << " - iterations : " << cnt;
if (cnt == 50)
{
ERROR << "Too much linesearch iterations ! We stop";
return ParabolicLinesearchReturnCode::MaximumIterations;
}
lambda_out = lambda;
switch (retcode) {
case 0:
return ParabolicLinesearchReturnCode::Success;
break;
case 1:
return ParabolicLinesearchReturnCode::LambdaTooSmall;
default:
return ParabolicLinesearchReturnCode::NotSupposedToHappen;
break;
}
}
template <class Program>
void ParabolicDriver<Program>::set_velocity(Vector& velocity_vector)
{
m_velocity.resize(program().get_tot_ndf());
m_velocity.setZero();
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
if (program().id_equation(dof) == no_equation)
{
m_velocity(dof) = velocity_vector(dof);
}
}
}
template <class Program>
void ParabolicDriver<Program>::reset_velocity()
{
for (index_t dof=0; dof<program().get_tot_ndf(); ++dof)
{
if (program().id_equation(dof) == no_equation) continue;
m_velocity(dof) = 0.0;
}
}
} // end namespace dfpmsolver
} // end namespace specmicp
c4science · Help