Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/kato.cpp
76105868d1eafeatures/netcdf-mpi

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>
#include <iostream>
/* -------------------------------------------------------------------------- */
__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>();
}
engine.registerNeumann();
}
/* -------------------------------------------------------------------------- */
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::reduce<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;
}
},
range<VectorProxy<Real, comp>>(*gap),
range<VectorProxy<Real, comp>>(*pressure));
g_static /= count_static != 0 ? count_static : 1;
g_static(comp - 1) = 0;
Loop::loop(
[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;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*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::loop(
[] CUDA_LAMBDA(Real& s, VectorProxy<Real, comp> sc) {
sc(comp - 1) = s;
},
surface, range<VectorProxy<Real, comp>>(*surfaceComp));
}
/* -------------------------------------------------------------------------- */
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);
// }
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureTresca() {
Loop::loop(
[this] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
if (p_N < 0) {
p_T = 0;
p(comp - 1) = 0;
} else if (p_T_norm > mu) { // TODO: replace name of variable mu
p_T *= mu / p_T_norm;
// p(comp - 1) += p_T_norm - mu;
}
},
range<VectorProxy<Real, comp>>(*pressure));
}
template void Kato::enforcePressureTresca<2>();
template void Kato::enforcePressureTresca<3>();
/* -------------------------------------------------------------------------- */
/**
* 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(bool use_tresca) {
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:
if (!use_tresca) {
beta = computeBeta<model_type::surface_1d>();
computeValuesForCost<model_type::surface_1d>(beta, lambda, eta, p_N, p_C);
} else {
computeValuesForCostTresca<model_type::surface_1d>(lambda, eta, p_N, p_C);
}
break;
case model_type::surface_2d:
if (!use_tresca) {
beta = computeBeta<model_type::surface_2d>();
computeValuesForCost<model_type::surface_2d>(beta, lambda, eta, p_N, p_C);
} else {
computeValuesForCostTresca<model_type::surface_2d>(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::reduce<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;
},
range<VectorProxy<Real, comp>>(*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::loop(
[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;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*gap), lambda, eta, p_N, p_C);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Kato::computeValuesForCostTresca(GridBase<Real>& lambda, GridBase<Real>& eta,
GridBase<Real>& p_N, GridBase<Real>& p_C) {
constexpr UInt comp = model_type_traits<type>::components;
Real g_N_min = Loop::reduce<operation::min>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
return g(comp - 1);
},
range<VectorProxy<Real, comp>>(*gap));
Real p_C_min = Loop::reduce<operation::min>(
[this] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
return mu - p_T_norm;
},
range<VectorProxy<Real, comp>>(*pressure));
Loop::loop(
[this, g_N_min, p_C_min] 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 - g_N_min;
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_N;
p_C_ = (p_N > 0) ? mu - p_T_norm - p_C_min : 0;
},
range<VectorProxy<Real, comp>>(*pressure),
range<VectorProxy<Real, comp>>(*gap), lambda, eta, p_N, p_C);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::computeFinalGap() {
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
Real g_N_min = Loop::reduce<operation::min>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
return g(comp - 1);
},
range<VectorProxy<Real, comp>>(*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