Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/polonsky_keer_tan.cpp
e8dcc4e0dab8features/netcdf-mpi

Last Change
Raw File

polonsky_keer_tan.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 "polonsky_keer_tan.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_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;
}
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;
*search_direction = 0.0;
Vector<Real, comp> gap_mean;
do {
enforcePressureMean<comp>(p0);
// Compute functional gradient
engine.solveNeumann(*pressure, *gap);
*gap -= surfaceComp;
computeMean<comp>(*gap, gap_mean, true);
*gap -= gap_mean;
// Compute search direction
Real R = computeSquaredNorm(*gap);
*search_direction *= 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>();
// *search_direction_backup = *gap;
Real alpha = tau;
while(enforcePressureConstraints<comp>(p0, alpha) &&
(cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations) {
cost = computeCost();
printState(n-1, cost, cost);
}
// Kato::enforcePressureConstraints<comp>(p0, 50);
// std::cout << "I got out at n = " << n << std::endl;
cost = computeCost();
printState(n, cost, cost);
} while ((cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
bool PolonskyKeerTan::enforcePressureConstraints(
GridBase<Real>& p0, Real& alpha) {
// *search_direction_backup = *gap;
// engine.solveNeumann(*pressure, *gap);
// *gap -= surfaceComp;
// Vector<Real, comp> gap_mean;
// computeMean<comp>(*gap, gap_mean, true);
// *gap -= gap_mean;
UInt na_count = Loop::stridedReduce<operation::plus>(
[alpha] CUDA_LAMBDA(VectorProxy<Real, comp>&& p,
VectorProxy<Real, comp>&& q,
VectorProxy<Real, comp>&& g) {
if (p(comp - 1) == 0.0 && g(comp - 1) < 0.0) {
Vector<Real, comp> _q = q;
_q *= alpha;
p -= _q;
return 1;
} else {
return 0;
}
},
*pressure, *gap, *gap);
// std::cout << na_count << std::endl;
if (na_count == 0) return false;
// Enforce mean value
enforcePressureMean<comp>(p0);
// Compute constraint gradient
computeGradientConstraints<comp>();
// Compute search direction
*search_direction = *gap;
truncateSearchDirection<comp>(true);
// Compute step size
alpha = computeStepSize<comp>(true);
// Update pressure
*search_direction *= alpha;
*pressure -= *search_direction;
enforcePressureCoulomb<comp>();
return true;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void PolonskyKeerTan::enforcePressureMean(GridBase<Real>& p0) {
Vector<Real, comp> pressure_mean;
computeMean<comp>(*pressure, pressure_mean, false);
*pressure -= pressure_mean;
addUniform<comp>(*pressure, p0);
// *pressure /= pressure_mean; // what if == 0 ?
// VectorProxy<Real, comp> _p0(p0(0));
// *pressure *= _p0;
// 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);
// p(comp - 1) *= p0(comp - 1) / pressure_mean(comp - 1);
// },
// *pressure);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void PolonskyKeerTan::computeMean(GridBase<Real>& field,
Vector<Real, comp>& mean,
bool on_c) {
UInt count = 0;
mean = Loop::stridedReduce<operation::plus>(
[&] CUDA_LAMBDA(VectorProxy<Real, comp>&& f,
VectorProxy<Real, comp>&& p) -> Vector<Real, comp> {
if ((!on_c) || p(comp - 1) > 0.0) {
count ++;
return f;
} else {
return 0;
}
},
field, *pressure);
mean /= count;
}
/* -------------------------------------------------------------------------- */
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::stridedLoop(
[] CUDA_LAMBDA(VectorProxy<Real, comp>&& t,
VectorProxy<Real, comp>&& p) {
if (p(comp - 1) == 0.0)
t = 0.0;
},
*search_direction, *pressure);
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
Real PolonskyKeerTan::computeStepSize(bool on_c) {
Vector<Real, comp> r_mean;
engine.solveNeumann(*search_direction, *projected_search_direction);
computeMean<comp>(*projected_search_direction, r_mean, on_c);
*projected_search_direction -= r_mean;
Real num = Loop::stridedReduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp>&& q, VectorProxy<Real, comp>&& t) {
Real dot = 0.0;
for (UInt i = 0; i < comp; i++) dot += q(i) * t(i);
return dot;
},
*gap, *search_direction);
Real denum = Loop::stridedReduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp>&& r, VectorProxy<Real, comp>&& t) {
Real dot = 0.0;
for (UInt i = 0; i < comp; i++) dot += r(i) * t(i);
return dot;
},
*projected_search_direction, *search_direction);
return num / denum;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void PolonskyKeerTan::computeGradientConstraints() {
engine.solveNeumann(*pressure, *gap);
*gap -= surfaceComp;
// *search_direction_backup = *gap;
// Vector<Real, comp> gap_mean;
// computeMean<comp>(*search_direction_backup, gap_mean, true);
// *search_direction_backup -= gap_mean;
Loop::stridedLoop(
[this] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
VectorProxy<Real, comp>&& p) {
VectorProxy<Real, comp - 1> g_T(g(0));
VectorProxy<Real, 1> g_N = g(comp - 1);
Real g_T_norm = g_T.l2norm();
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
g_N += mu * g_T_norm;
if (p_T_norm == 0.0) {
g_T *= 0.0;//-1.0;
} else {
g_T = p_T;
g_T *= - g_T_norm / p_T_norm;
}
},
*gap, *pressure);
Vector<Real, comp> gap_mean;
computeMean<comp>(*gap, gap_mean, true);
*gap -= gap_mean;
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help