Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/solvers/kato.cpp
1c342c8abb54maxwell_viscoelastic…

Last Change
Raw File

kato.cpp
View Options

/**
* @file
*
* @author Son Pham-Ba <son.phamba@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "kato.hh"
#include "elastic_functional.hh"
#include "loop.hh"
#include "fft_plan_manager.hh"
#include <iomanip>
#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Kato::Kato(Model& model, const GridBase<Real>& surface, Real tolerance,
Real mu)
: ContactSolver(model, surface, tolerance),
engine(model.getBEEngine()), mu(mu) {
if (model.getType() != model_type::surface_1d &&
model.getType() != model_type::surface_2d) {
TAMAAS_EXCEPTION(
"Model type is not compatible with Kato solver");
}
gap = this->_gap.get(); // locally allocated
pressure = &model.getTraction();
N = pressure->getNbPoints();
if (model.getType() == model_type::surface_1d) {
initSurfaceWithComponents<model_type::surface_1d>();
} else {
initSurfaceWithComponents<model_type::surface_2d>();
}
}
/* -------------------------------------------------------------------------- */
Real Kato::solve(GridBase<Real>& p0, UInt proj_iter) {
if (p0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
"Target mean pressure does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveTmpl<model_type::surface_1d>(p0, proj_iter);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0, proj_iter);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveTmpl(GridBase<Real>& p0, UInt proj_iter) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
pressure->uniformSetComponents(p0);
do {
computeGradient<comp>();
*pressure -= *gap;
enforcePressureConstraints<comp>(p0, proj_iter);
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
Real Kato::solveRelaxed(GridBase<Real>& g0) {
if (g0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
"Target mean gap does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveRelaxedTmpl<model_type::surface_1d>(g0);
break;
case model_type::surface_2d:
cost = solveRelaxedTmpl<model_type::surface_2d>(g0);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveRelaxedTmpl(GridBase<Real>& g0) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
*pressure = 0;
do {
engine.solveNeumann(*pressure, *gap);
addUniform<comp>(*gap, g0);
*gap -= *surfaceComp;
*pressure -= *gap;
enforcePressureCoulomb<comp>();
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
Real Kato::solveRegularized(GridBase<Real>& p0, Real r) {
if (p0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
"Target mean pressure does not have the right number of components");
}
Real cost = 0;
switch (model.getType()) {
case model_type::surface_1d:
cost = solveRegularizedTmpl<model_type::surface_1d>(p0, r);
break;
case model_type::surface_2d:
cost = solveRegularizedTmpl<model_type::surface_2d>(p0, r);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Kato::solveRegularizedTmpl(GridBase<Real>& p0, Real r) {
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0;
UInt n = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
pressure->uniformSetComponents(p0);
do {
// enforcePressureMean<comp>(p0);
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
// Impose zero tangential displacement in non-sliding zone
UInt count_static = 0;
Vector<Real, comp> g_static = Loop::stridedReduce<operation::plus>(
[&] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
VectorProxy<Real, comp>&& p) -> Vector<Real, comp> {
VectorProxy<Real, comp - 1> p_T = p(0);
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact
count_static ++;
return g; // to compute mean of g_T
} else {
return 0;
}
},
*gap, *pressure);
g_static /= count_static != 0 ? count_static : 1;
g_static(comp - 1) = 0;
Loop::stridedLoop(
[this, r, g_static] CUDA_LAMBDA(VectorProxy<Real, comp>&& p,
VectorProxy<Real, comp>&& g) {
// Add frictional term to gradient of functional
g -= g_static;
Vector<Real, comp> _g = g; // copy
VectorProxy<Real, comp - 1> g_T (g(0));
VectorProxy<Real, 1> g_N (g(comp - 1));
Real g_T_norm = g_T.l2norm();
// g_N += mu * regularize(g_T_norm, r) * g_T_norm;
g_N += mu * g_T_norm;
// Update pressure with gradient
// _g *= 0.1;
p -= g;
// Truncate negative normal pressure
VectorProxy<Real, comp - 1> p_T (p(0));
VectorProxy<Real, 1> p_N (p(comp - 1));
if (p_N(0) < 0) p_N = 0;
// Set tangential pressure
p_T = g_T;
if (g_T_norm != 0)
p_T *= - mu * p_N(0) * regularize(g_T_norm, r) / g_T_norm;
},
*pressure, *gap);
// enforcePressureMean<comp>(p0);
// enforcePressureCoulomb<comp>();
enforcePressureConstraints<comp>(p0, 50);
cost = computeCost();
printState(n, cost, cost);
} while (std::abs(cost) > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Kato::initSurfaceWithComponents() {
constexpr UInt comp = model_type_traits<type>::components;
surfaceComp =
allocateGrid<true, Real>(type, model.getDiscretization(), comp);
*surfaceComp = 0;
Loop::stridedLoop(
[] CUDA_LAMBDA(Real& s, VectorProxy<Real, comp>&& sc) {
sc(comp - 1) = s;
},
surface, *surfaceComp);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::computeGradient() {
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
// Impose zero tangential displacement in non-sliding zone
UInt count_static = 0;
Vector<Real, comp> g_static = Loop::stridedReduce<operation::plus>(
[&] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
VectorProxy<Real, comp>&& p) -> Vector<Real, comp> {
VectorProxy<Real, comp - 1> p_T = p(0);
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact
count_static ++;
return g; // to compute mean of g_T
} else {
return 0;
}
},
*gap, *pressure);
if (count_static != 0) {
g_static /= count_static;
// g_static(comp - 1) = 0;
} else { // if no sticking zone, mean computed on sliding zone
count_static = 0;
Real g_N_static = Loop::stridedReduce<operation::plus>(
[&] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
VectorProxy<Real, comp>&& p) {
if (p(comp - 1) > 0) {
count_static ++;
return g(comp - 1);
} else {
return 0.0;
}
},
*gap, *pressure);
g_static(comp - 1) = g_N_static / count_static;
}
// Add frictionnal term to functional
Loop::stridedLoop(
[this, g_static] CUDA_LAMBDA(VectorProxy<Real, comp>&& g) {
g -= g_static;
VectorProxy<Real, comp - 1> g_T = g(0);
Real g_T_norm = g_T.l2norm();
g(comp - 1) += mu * g_T_norm; // Frictionnal work
},
*gap);
// Loop::stridedLoop(
// [this] CUDA_LAMBDA(VectorProxy<Real, comp>&& g, VectorProxy<Real, comp>&& p) {
// VectorProxy<Real, comp - 1> g_T = g(0);
// Real g_T_norm = g_T.l2norm();
// VectorProxy<Real, 1> g_N = g(comp - 1);
// VectorProxy<Real, comp - 1> p_T = p(0);
// Real p_T_norm = p_T.l2norm();
// VectorProxy<Real, 1> p_N = p(comp - 1);
// if (p_T_norm != 0) {
// g_T = p_T;
// g_T *= -g_T_norm / p_T_norm;
// }
// g_N += mu * g_T_norm; // Frictionnal work
// },
// *gap, *pressure);
}
/* -------------------------------------------------------------------------- */
/**
* Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$.
*/
template <UInt comp>
void Kato::enforcePressureConstraints(GridBase<Real>& p0, UInt proj_iter) {
for (UInt i = 0; i < proj_iter; i++) {
enforcePressureMean<comp>(p0);
enforcePressureCoulomb<comp>();
}
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureMean(GridBase<Real>& p0) {
Vector<Real, comp> corr = computeMean<comp>(*pressure);
VectorProxy<Real, comp> _p0 = p0(0);
corr -= _p0;
*pressure -= corr;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureCoulomb() {
Loop::stridedLoop(
[this] CUDA_LAMBDA(VectorProxy<Real, comp>&& p) {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_sqrd= p_T.l2squared();
// Projection normale au cône de friction
bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N);
bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu);
if (cond2) {
p_T = 0;
p(comp - 1) = 0;
} else if (!cond1) {
Real p_T_norm = std::sqrt(p_T_sqrd);
Real k = (p_N + mu * p_T_norm) / (1 + mu * mu);
p_T *= k * mu / p_T_norm;
p(comp - 1) = k;
}
},
*pressure);
}
/* -------------------------------------------------------------------------- */
/**
* Compute mean of the field taking each component separately.
*/
template <UInt comp>
Vector<Real, comp> Kato::computeMean(GridBase<Real>& field) {
Vector<Real, comp> mean = Loop::stridedReduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp>&& f) -> Vector<Real, comp> {
return f;
},
field);
mean /= N;
return mean;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::addUniform(GridBase<Real>& field, GridBase<Real>& vec) {
VectorProxy<Real, comp> _vec(vec(0));
field += _vec;
}
/* -------------------------------------------------------------------------- */
Real Kato::computeCost() {
UInt N = pressure->getNbPoints();
Real beta = 0;
Grid<Real, 1> lambda({N}, 1);
Grid<Real, 1> eta({N}, 1);
Grid<Real, 1> p_N({N}, 1);
Grid<Real, 1> p_C({N}, 1);
switch (model.getType()) {
case model_type::surface_1d:
beta = computeBeta<model_type::surface_1d>();
computeValuesForCost<model_type::surface_1d>(beta, lambda, eta, p_N, p_C);
break;
case model_type::surface_2d:
beta = computeBeta<model_type::surface_2d>();
computeValuesForCost<model_type::surface_2d>(beta, lambda, eta, p_N, p_C);
break;
default:
break;
}
return p_N.dot(lambda) + p_C.dot(eta);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
Real Kato::computeBeta() {
constexpr UInt comp = model_type_traits<type>::components;
return Loop::stridedReduce<operation::max>(
[this] CUDA_LAMBDA(VectorProxy<Real, comp>&& g) {
VectorProxy<Real, comp - 1> g_T(g(0));
Real g_N = g(comp - 1);
Real g_T_norm = g_T.l2norm();
return mu * g_T_norm - g_N;
},
*gap);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Kato::computeValuesForCost(Real beta, GridBase<Real>& lambda, GridBase<Real>& eta,
GridBase<Real>& p_N, GridBase<Real>& p_C) {
constexpr UInt comp = model_type_traits<type>::components;
Loop::stridedLoop(
[this, beta] CUDA_LAMBDA(VectorProxy<Real, comp>&& p,
VectorProxy<Real, comp>&& g,
Real& lambda_,
Real& eta_,
Real& p_N_,
Real& p_C_) {
VectorProxy<Real, comp - 1> g_T(g(0));
Real g_N = g(comp - 1);
Real g_T_norm = g_T.l2norm();
lambda_ = g_N - mu * g_T_norm + beta;
eta_ = g_T_norm;
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
p_N_ = p(comp - 1);
p_C_ = mu * p_N - p_T_norm;
},
*pressure, *gap, lambda, eta, p_N, p_C);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::computeFinalGap() {
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
Real g_N_min = Loop::stridedReduce<operation::min>(
[] CUDA_LAMBDA(VectorProxy<Real, comp>&& g) {
return g(comp - 1);
},
*gap);
Grid<Real, 1> g_shift({comp}, 1);
g_shift = 0;
g_shift(comp - 1) = -g_N_min;
*gap += *surfaceComp;
addUniform<comp>(*gap, g_shift);
model.getDisplacement() = *gap;
}
/* -------------------------------------------------------------------------- */
Real Kato::regularize(Real x, Real r) {
Real xr = x / r;
return xr / (1 + std::abs(xr));
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help