Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/epic.cpp
2316c9ebcdd7features/netcdf-mpi

Last Change
Raw File

epic.cpp
View Options

/**
* @file
* @section 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 "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());
#define SET_VIEWS(_, __, type) \
case type: { \
setViews<type>(); \
break; \
}
switch (model.getType()) {
BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
}
#undef SET_VIEWS
}
/* -------------------------------------------------------------------------- */
void EPICSolver::solve(std::vector<Real> load) {
UInt n = 0, nmax = 1000;
Real error = 0.;
const GridBase<Real> initial_surface(surface);
Real normalizing_factor = std::sqrt(initial_surface.var());
GridBase<Real> previous_residual(*residual_disp);
GridBase<Real> relaxation_direction(*residual_disp);
previous_residual = 0;
pressure = *pressure_inc; // set previous pressure to current model pressure
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 << std::endl;
} 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();
}
/* -------------------------------------------------------------------------- */
void EPICSolver::acceleratedSolve(std::vector<Real> load) {
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 << std::endl;
} 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();
}
/* -------------------------------------------------------------------------- */
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(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
c4science · Help