Page MenuHomec4science
Files F122571917

epic.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 20, 13:33

epic.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 "epic.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
EPICSolver::EPICSolver(ContactSolver& csolver, EPSolver& epsolver,
Real tolerance, Real relaxation)
: csolver(csolver), epsolver(epsolver), tolerance(tolerance),
relaxation(relaxation) {
auto& model = csolver.getModel();
surface.wrap(csolver.getSurface());
pressure.resize(model.getTraction().getNbPoints());
model_type_dispatch(
[this](auto&& _) {
constexpr auto type = std::decay_t<decltype(_)>::value;
setViews<type>();
},
model.getType());
}
/* -------------------------------------------------------------------------- */
Real EPICSolver::solve(const std::vector<Real>& load) {
const GridBase<Real> initial_surface{surface};
GridBase<Real> previous_residual{*residual_disp},
relaxation_direction{*residual_disp};
UInt n = 0;
Real error = 0, normalizing_factor = std::sqrt(initial_surface.var());
pressure = *pressure_inc; // set previous pressure to current model pressure
previous_residual = 0; // initial elastic guess
do {
fixedPoint(relaxation_direction, previous_residual, initial_surface, load);
relaxation_direction -= previous_residual;
relaxation_direction *= relaxation;
error = computeError(*residual_disp, previous_residual, normalizing_factor);
previous_residual += relaxation_direction;
Logger().get(LogLevel::info) << "[EPIC] 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;
}
/* -------------------------------------------------------------------------- */
Real EPICSolver::acceleratedSolve(const std::vector<Real>& load) {
Logger().get(LogLevel::warning)
<< "Irons & Truck iteration often has convergence issues. Using the "
"AndersonMixing solver is encouraged\n";
UInt n = 0, nmax = 1000;
Real error = 0.;
const GridBase<Real> initial_surface(surface);
Real normalizing_factor = std::sqrt(initial_surface.var());
GridBase<Real> doubleg(*residual_disp), dg(*residual_disp),
d2x(*residual_disp), x(*residual_disp), x_prev(*residual_disp);
do {
// First compute g(x) and g(g(x))
fixedPoint(dg, x, initial_surface, load);
fixedPoint(doubleg, dg, initial_surface, load);
// dg contains g(x): compute delta x
d2x = dg;
d2x -= x;
// compute delta g(x)
dg *= -1;
dg += doubleg;
// now we have delta g(x) and delta x, compute delta^2 x
d2x *= -1;
d2x += dg;
// Apply Irons & Truck iteration
auto f = dg.dot(d2x) / d2x.dot(d2x);
dg *= -f;
x = doubleg;
x += dg;
error = computeError(x, x_prev, normalizing_factor);
x_prev = x;
Logger().get(LogLevel::info) << "[EPIC] error = " << std::scientific
<< error << std::fixed << "\n";
} while (error > tolerance && n++ < nmax);
// computing full pressure
pressure += *pressure_inc;
// setting the model pressure to full converged pressure
*pressure_inc = pressure;
// updating plastic state
epsolver.updateState();
return error;
}
/* -------------------------------------------------------------------------- */
void EPICSolver::fixedPoint(GridBase<Real>& result, const GridBase<Real>& x,
const GridBase<Real>& initial_surface,
std::vector<Real> load) {
surface = initial_surface;
surface -= x;
csolver.solve(std::move(load));
*pressure_inc -= pressure;
epsolver.solve();
result = *residual_disp;
}
/* -------------------------------------------------------------------------- */
Real EPICSolver::computeError(const GridBase<Real>& current,
const GridBase<Real>& prev, Real factor) const {
GridBase<Real> error(current);
error -= prev;
return std::sqrt(error.dot(error)) / factor;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas

Event Timeline

c4science · Help