Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/solvers/polonsky_keer_tan.cpp
531937f6e06ffeatures/declarative

Last Change
Raw File

polonsky_keer_tan.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 "polonsky_keer_tan.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
namespace tamaas {
PolonskyKeerTan::PolonskyKeerTan(Model& model, const GridBase<Real>& surface,
Real tolerance, Real mu)
: Kato(model, surface, tolerance, mu) {
search_direction =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
search_direction_backup =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
projected_search_direction =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerTan::solve(GridBase<Real>& p0) {
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);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0);
break;
default:
break;
}
return cost;
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerTan::solveTresca(GridBase<Real>& p0) {
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, true);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0, true);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real PolonskyKeerTan::solveTmpl(GridBase<Real>& p0, bool use_tresca) {
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
constexpr UInt comp = model_type_traits<type>::components;
Real cost = 0.0;
Real error = 0.0;
UInt n = 0;
pressure->uniformSetComponents(p0);
Real R_old = 1.0;
*search_direction = 0.0;
do {
// Enforce external condition (should be at the end)
enforcePressureMean<comp>(p0);
// Compute functional gradient
computeGradient<comp>(use_tresca);
// Compute search direction
Real R = computeSquaredNorm(*gap);
*search_direction *= R / R_old;
*search_direction += *gap;
R_old = R;
// Compute step size
Real tau = computeStepSize<comp>(false);
// Update pressure
*search_direction *= tau;
*pressure -= *search_direction;
*search_direction *= model.getSystemSize()[0] / model.getYoungModulus();
// *search_direction *= 1/tau; // this line should *theoreticaly* replace the one above
// Enforce constraints
if (!use_tresca) {
enforcePressureCoulomb<comp>();
} else {
enforcePressureTresca<comp>();
}
cost = computeCost(use_tresca);
error = cost / pressure->getNbPoints() / model.getSystemSize()[0] / model.getYoungModulus();
printState(n, cost, error);
} while (error > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return error;
}
// Original algorithm
// template <model_type type>
// Real PolonskyKeerTan::solveTmpl(GridBase<Real>& p0) {
// // Printing column headers
// std::cout << std::setw(5) << "Iter"
// << " " << std::setw(15) << "Cost_f"
// << " " << std::setw(15) << "Error" << '\n'
// << std::fixed;
// constexpr UInt comp = model_type_traits<type>::components;
// Real cost = 0.0;
// UInt n = 0;
// pressure->uniformSetComponents(p0);
// Real R_old = 1.0;
// Real delta = 0.0;
// *search_direction = 0.0;
// do {
// // Enforce external condition (should be at the end)
// enforcePressureMean<comp>(p0);
// // Compute functional gradient
// computeGradient<comp>();
// Vector<Real, comp> gap_mean = computeMean<comp>(*gap, true);
// for (UInt i = 0; i < comp - 1; i++) gap_mean(i) = 0;
// *gap -= gap_mean;
// // Compute search direction
// Real R = computeSquaredNorm(*gap);
// *search_direction *= delta * R / R_old;
// *search_direction += *gap;
// R_old = R;
// truncateSearchDirection<comp>(true);
// // Compute step size
// Real tau = computeStepSize<comp>(true);
// // Update pressure
// *search_direction *= tau;
// *pressure -= *search_direction;
// // Enforce constraints
// enforcePressureCoulomb<comp>();
// // Empty set of inadmissible gaps
// UInt na_count = Loop::stridedReduce<operation::plus>(
// [tau] CUDA_LAMBDA(VectorProxy<Real, comp>&& p,
// VectorProxy<Real, comp>&& g) {
// if (p(comp - 1) == 0.0 && g(comp - 1) < 0.0) {
// Vector<Real, comp> _g = g;
// _g *= tau;
// p -= _g;
// return 1;
// } else {
// return 0;
// }
// },
// *pressure, *gap);
// delta = (na_count > 0) ? 0.0 : 1.0;
// // Enforce external condition
// // enforcePressureMean<comp>(p0);
// cost = computeCost();
// printState(n, cost, cost);
// } while (cost > this->tolerance && n++ < this->max_iterations);
// computeFinalGap<comp>();
// return cost;
// }
/* -------------------------------------------------------------------------- */
template <UInt comp>
void PolonskyKeerTan::enforcePressureMean(GridBase<Real>& p0) {
Vector<Real, comp> pressure_mean =
computeMean<comp>(*pressure, false);
// *pressure -= pressure_mean;
// addUniform<comp>(*pressure, p0);
// for (UInt i = 0; i < comp; i++)
// if (pressure_mean(i) == 0) pressure_mean(i) = 1.0;
// *pressure /= pressure_mean;
// VectorProxy<Real, comp> _p0(p0(0));
// *pressure *= _p0;
Loop::loop(
[pressure_mean, p0] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
for (UInt i = 0; i < comp - 1; i++)
p(i) += p0(i) - pressure_mean(i);
p(comp - 1) *= p0(comp - 1) / pressure_mean(comp - 1);
},
range<VectorProxy<Real, comp>>(*pressure));
// Loop::stridedLoop(
// [pressure_mean, p0] CUDA_LAMBDA(VectorProxy<Real, comp>&& p) {
// for (UInt i = 0; i < comp - 1; i++)
// p(i) *= p0(i) / (pressure_mean(i) != 0 ? pressure_mean(i) : 1);
// p(comp - 1) += p0(comp - 1) - pressure_mean(comp - 1);
// },
// *pressure);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
Vector<Real, comp> PolonskyKeerTan::computeMean(GridBase<Real>& field,
bool on_c) {
UInt count = Loop::reduce<operation::plus>(
[on_c] CUDA_LAMBDA(VectorProxy<Real, comp> p) -> UInt {
if ((!on_c) || p(comp - 1) > 0.0) {
return 1;
} else {
return 0;
}
},
range<VectorProxy<Real, comp>>(*pressure));
Vector<Real, comp> mean = Loop::reduce<operation::plus>(
[on_c] CUDA_LAMBDA(VectorProxy<Real, comp> f,
VectorProxy<Real, comp> p) -> Vector<Real, comp> {
if ((!on_c) || p(comp - 1) > 0.0) {
return f;
} else {
return 0;
}
},
range<VectorProxy<Real, comp>>(field), range<VectorProxy<Real, comp>>(*pressure));
mean /= count;
return mean;
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerTan::computeSquaredNorm(GridBase<Real>& field) {
return Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real& f) {
return f * f;
},
field);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void PolonskyKeerTan::truncateSearchDirection(bool on_c) {
if (!on_c) return;
Loop::loop(
[] CUDA_LAMBDA(VectorProxy<Real, comp> t, VectorProxy<Real, comp> p) {
if (p(comp - 1) == 0.0)
t = 0.0;
},
range<VectorProxy<Real, comp>>(*search_direction),
range<VectorProxy<Real, comp>>(*pressure));
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
Real PolonskyKeerTan::computeStepSize(bool on_c) {
engine.solveNeumann(*search_direction, *projected_search_direction);
Vector<Real, comp> r_mean =
computeMean<comp>(*projected_search_direction, on_c);
*projected_search_direction -= r_mean;
Real num = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> q,
VectorProxy<Real, comp> t) {
return q.dot(t);
},
range<VectorProxy<Real, comp>>(*gap),
range<VectorProxy<Real, comp>>(*search_direction));
Real denum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> r,
VectorProxy<Real, comp> t) {
return r.dot(t);
},
range<VectorProxy<Real, comp>>(*projected_search_direction),
range<VectorProxy<Real, comp>>(*search_direction));
// if (denum == 0) denum = 1.0;
return num / denum;
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
c4science · Help