diff --git a/src/solvers/anderson.cpp b/src/solvers/anderson.cpp
index 529d1eb..68739c0 100644
--- a/src/solvers/anderson.cpp
+++ b/src/solvers/anderson.cpp
@@ -1,192 +1,190 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "anderson.hh"
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
AndersonMixing::AndersonMixing(ContactSolver& csolver, EPSolver& epsolver,
Real tolerance, UInt memory)
: EPICSolver(csolver, epsolver, tolerance, 0.5), M(memory) {}
/* -------------------------------------------------------------------------- */
Real AndersonMixing::solve(const std::vector& load) {
UInt n = 0, nmax = 1000;
Real error = 0.;
const GridBase initial_surface(surface);
Real normalizing_factor = std::sqrt(initial_surface.var());
- GridBase x(*residual_disp), x_prev(*residual_disp);
+ GridBase x{*residual_disp}, x_prev{*residual_disp};
memory_t memory, residual_memory;
do {
{
// Compute a residual vector for this iteration
- GridBase F(*residual_disp);
+ GridBase 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
if (memory.size() > M + 1)
memory.pop_back();
if (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 << std::endl;
} while (error > tolerance && n++ < nmax);
- std::cout << "Converged in " << n << " iterations\n";
-
// 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 AndersonMixing::mixingUpdate(GridBase x,
std::vector 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 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{A.sizes(), 1}, Grid{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> LU, Grid 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 AndersonMixing::computeGamma(const memory_t& residual) {
- auto N = residual.size() - 1;
- Grid A({N, N}, 1);
- Grid b({N}, 1);
+ UInt N = residual.size() - 1;
+ Grid A{{N, N}, 1};
+ Grid 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(A), b);
+ auto x = LUsubstitute(factorize(std::move(A)), std::move(b));
return std::vector(x.begin(), x.end());
}
} // namespace tamaas