Page MenuHomec4science
Files F92106768

anderson.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Nov 17, 10:08

anderson.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2024 Lucas Frérot
*
* 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 "anderson.hh"
#include <algorithm>
#include <iomanip>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
AndersonMixing::AndersonMixing(ContactSolver& csolver, EPSolver& epsolver,
Real tolerance, UInt memory)
: EPICSolver(csolver, epsolver, tolerance, 1), M(memory) {}
/* -------------------------------------------------------------------------- */
Real AndersonMixing::solve(const std::vector<Real>& load) {
const GridBase<Real> initial_surface{surface};
GridBase<Real> x{*residual_disp}, x_prev{*residual_disp};
memory_t memory, residual_memory;
UInt n = 0;
Real error = 0, normalizing_factor = std::sqrt(initial_surface.var());
pressure = *pressure_inc; // set previous pressure to current model pressure
x = x_prev = 0; // initial elastic guess
do {
{
// Compute a residual vector for this iteration
GridBase<Real> F{*residual_disp};
fixedPoint(F, x, initial_surface, load);
F -= x; // F = g(x) - x
memory.push_front(x);
residual_memory.push_front(std::move(F));
}
// Keep M previous iterations: memory has M + 1 vectors
while (memory.size() > M + 1)
memory.pop_back();
while (residual_memory.size() > M + 1)
residual_memory.pop_back();
x = mixingUpdate(std::move(x), computeGamma(residual_memory), memory,
residual_memory, relaxation);
error = computeError(x, x_prev, normalizing_factor);
x_prev = x;
Logger().get(LogLevel::info) << "[Anderson] error = " << std::scientific
<< error << std::fixed << "\n";
} while (error > tolerance && n++ < max_iterations);
// computing full pressure
pressure += *pressure_inc;
// setting the model pressure to full converged pressure
*pressure_inc = pressure;
// updating plastic state
epsolver.updateState();
return error;
}
/* -------------------------------------------------------------------------- */
GridBase<Real> AndersonMixing::mixingUpdate(GridBase<Real> x,
std::vector<Real> gamma,
const memory_t& memory,
const memory_t& residual_memory,
Real relaxation) {
Loop::loop([beta = relaxation](Real& x, Real Fl) { x += beta * Fl; }, x,
residual_memory.front());
for (UInt m = 1; m < memory.size(); ++m) {
Loop::loop(
[gm = gamma[m - 1], beta = relaxation](Real& x, Real xmp1, Real xm,
Real Fmp1, Real Fm) {
// Eyert (10.1006/jcph.1996.0059) p282 eq (7.7)
x -= gm * ((xmp1 - xm) + beta * (Fmp1 - Fm));
},
x, memory[m - 1], memory[m], residual_memory[m - 1],
residual_memory[m]);
}
return x;
}
/* -------------------------------------------------------------------------- */
/// Crout(e) (Antoine... et Daniel...) factorization from
/// https://en.wikipedia.org/wiki/Crout_matrix_decomposition
auto factorize(Grid<Real, 2> A) {
TAMAAS_ASSERT(A.sizes()[0] == A.sizes()[1], "Matrix is not square");
const auto N = A.sizes()[0];
auto LU =
std::make_pair(Grid<Real, 2>{A.sizes(), 1}, Grid<Real, 2>{A.sizes(), 1});
auto& L = LU.first;
auto& U = LU.second;
for (UInt i = 0; i < N; ++i) {
U(i, i) = 1;
}
for (UInt j = 0; j < N; ++j) {
for (UInt i = j; i < N; ++i) {
double sum = 0;
for (UInt k = 0; k < j; ++k) {
sum += L(i, k) * U(k, j);
}
L(i, j) = A(i, j) - sum;
}
for (UInt i = j; i < N; ++i) {
double sum = 0;
for (UInt k = 0; k < j; ++k) {
sum += L(j, k) * U(k, i);
}
U(j, i) = (A(j, i) - sum) / L(j, j);
}
}
return LU;
}
/* -------------------------------------------------------------------------- */
auto LUsubstitute(std::pair<Grid<Real, 2>, Grid<Real, 2>> LU, Grid<Real, 1> b) {
const auto& L = LU.first;
const auto& U = LU.second;
const auto N = b.sizes()[0];
// Forward substitution
for (UInt m = 0; m < N; ++m) {
for (UInt i = 0; i < m; ++i)
b(m) -= L(m, i) * b(i);
b(m) /= L(m, m);
}
// Backward substitution
for (UInt l = 0; l < N; ++l) {
auto m = N - 1 - l;
for (UInt i = m + 1; i < N; ++i)
b(m) -= U(m, i) * b(i);
b(m) /= U(m, m);
}
return b;
}
/* -------------------------------------------------------------------------- */
std::vector<Real> AndersonMixing::computeGamma(const memory_t& residual) {
UInt N = residual.size() - 1;
Grid<Real, 2> A{{N, N}, 1};
Grid<Real, 1> b{{N}, 1};
const auto& Fl = residual.front();
// Fill RHS vector
for (UInt n = 1; n < residual.size(); ++n) {
b(n - 1) = residual[n - 1].dot(Fl) - residual[n].dot(Fl);
}
/// Regularization from Eyert
auto w0 = 0;
// Fill matrix
for (UInt n = 1; n < residual.size(); ++n) {
for (UInt m = 1; m < residual.size(); ++m) {
A(n - 1, m - 1) = residual[n - 1].dot(residual[m - 1]) -
residual[n - 1].dot(residual[m]) -
residual[n].dot(residual[m - 1]) +
residual[n].dot(residual[m]);
A(n - 1, m - 1) *= (1 + w0 * w0 * (n == m));
}
}
auto x = LUsubstitute(factorize(std::move(A)), std::move(b));
return std::vector<Real>(x.begin(), x.end());
}
} // namespace tamaas

Event Timeline

c4science · Help