Homec4science
Diffusion tamaas (petsc) Browse

tamaas/src/solvers/dfsane_solver.cpp
9fe244aa2311petsc

Last Change
Raw File

dfsane_solver.cpp
View Options

/**
* @file
* LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "dfsane_solver.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
DFSANESolver::DFSANESolver(Residual& residual)
: EPSolver(residual),
search_direction(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
previous_residual(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
current_x(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
delta_x(residual.getVector().dataSize(),
residual.getVector().getNbComponents()),
delta_residual(residual.getVector().dataSize(),
residual.getVector().getNbComponents()) {}
void DFSANESolver::solve() {
EPSolver::beforeSolve();
auto& x_prev = *_x;
UInt k = 0, nmax = 500;
Real F_norm = 0;
previous_merits.clear();
_residual.computeResidual(x_prev);
previous_residual = _residual.getVector();
const auto F_norm0 = previous_residual.l2norm();
eta = [F_norm0](UInt k) { return F_norm0 / ((1 + k) * (1 + k)); };
const std::pair<Real, Real> sigma_extrema{1e-10, 1e10};
Real sigma = 1;
do {
computeSearchDirection(sigma);
lineSearch(eta(k));
// No need for a solution update
// it was done automatically in line search
sigma = computeSpectralCoeff(sigma_extrema);
x_prev = current_x;
previous_residual = _residual.getVector();
F_norm = previous_residual.l2norm();
Logger().get(LogLevel::debug) << TAMAAS_DEBUG_MSG(k << ' ' << F_norm);
} while (F_norm > getTolerance() and k++ < nmax);
if (k >= nmax)
TAMAAS_EXCEPTION("DF-SANE did not converge");
Logger().get(LogLevel::info) << "DF-SANE: sucessful convergence (" << k
<< " iterations, " << getTolerance() << ")\n";
_residual.computeResidualDisplacement(*_x);
}
Real DFSANESolver::computeSpectralCoeff(const std::pair<Real, Real>& bounds) {
delta_x = current_x;
delta_x -= *_x;
delta_residual = _residual.getVector();
delta_residual -= previous_residual;
auto trial = delta_x.dot(delta_x) / delta_x.dot(delta_residual);
if (std::abs(trial) >= bounds.first and std::abs(trial) <= bounds.second)
return trial;
else {
auto F = previous_residual.l2norm();
// Rule from Cruz et al. (2006)
if (F > 1)
return 1;
else if (F < 1e-5)
return 1e5;
else
return 1 / F;
}
}
void DFSANESolver::computeSearchDirection(Real sigma) {
_residual.computeResidual(*_x);
search_direction = _residual.getVector();
search_direction *= -sigma;
}
void DFSANESolver::lineSearch(Real eta_k) {
const UInt M = 10;
const Real gamma = 1e-4;
Real alphap = 1, alpham = 1;
const std::pair<Real, Real> tau_extrema{0.1, 0.5};
// Merit function (nexp = 2)
const auto f = [&](Real alpha, Real sign) {
current_x = search_direction;
current_x *= alpha * sign;
current_x += *_x;
_residual.computeResidual(current_x);
return _residual.getVector().dot(_residual.getVector());
};
// Manage previous merits
const auto fk = _residual.getVector().dot(_residual.getVector());
previous_merits.push_front(fk);
if (previous_merits.size() > M)
previous_merits.pop_back();
const auto fbar =
*std::max_element(previous_merits.begin(), previous_merits.end());
do {
Real fp = f(alphap, +1);
if (fp < fbar + eta_k - gamma * alphap * alphap * fk)
return;
Real fm = f(alpham, -1);
if (fm < fbar + eta_k - gamma * alpham * alpham * fk)
return;
alphap = computeAlpha(alphap, fp, fk, tau_extrema);
alpham = computeAlpha(alpham, fm, fk, tau_extrema);
if (std::isnan(alphap) or std::isnan(alpham))
TAMAAS_EXCEPTION("NaN error in line search");
} while (true);
}
Real DFSANESolver::computeAlpha(const Real alpha, Real f, Real fk,
const std::pair<Real, Real>& bounds) {
const auto alphat = alpha * alpha * fk / (f + (2 * alpha - 1) * fk);
if (alphat < bounds.first * alpha)
return bounds.first * alpha;
else if (alphat > bounds.second * alpha)
return bounds.second * alpha;
else
return alphat;
}
} // namespace tamaas
c4science · Help