Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/solvers/polonsky_keer_rey.cpp
b6fe91b593c4archive/legacy_mpi

Last Change
Raw File

polonsky_keer_rey.cpp
View Options

/**
* @file
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-2017 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_rey.hh"
#include "elastic_functional.hh"
#include "loop.hh"
#include "model_type.hh"
#include <iomanip>
#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
PolonskyKeerRey::PolonskyKeerRey(Model& model, GridBase<Real>& surface,
Real tolerance, type variable_type,
type constraint_type)
: ContactSolver(model, surface, tolerance), variable_type(variable_type),
constraint_type(constraint_type), functional(model.getBEEngine()) {
if (model.getType() != model_type::basic_1d &&
model.getType() != model_type::basic_2d) {
TAMAAS_EXCEPTION(
"Model type is not compatible with PolonskyKeerRey solver");
}
/// Gap is locally allocated
auto gap_ =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
auto pressure_ = &model.getTraction();
switch (variable_type) {
case gap:
primal = gap_;
dual = pressure_;
functional.addFunctionalTerm(
std::make_unique<functional::ElasticFunctionalGap>(model.getBEEngine(),
surface));
break;
case pressure:
primal = pressure_;
dual = gap_;
functional.addFunctionalTerm(
std::make_unique<functional::ElasticFunctionalPressure>(
model.getBEEngine(), surface));
break;
}
search_direction =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
projected_search_direction =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerRey::solve(Real target) {
Real G = 0, Gold = 1, error = 0, error_norm = 1;
UInt n = 0;
bool delta = false;
*search_direction = 0;
//*dual = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
// Setting uniform value if constraint
if (constraint_type == variable_type && std::abs(primal->sum()) <= 0)
*primal = target;
else if (constraint_type != variable_type)
*primal = this->surface_stddev;
do {
// Computing functional gradient
functional.computeGradF(*primal, *dual);
// Enforcing dual constraint via gradient
if (constraint_type != variable_type) {
// Real dbar = meanOnUnsaturated(*dual);
// *dual += 2*target + dbar;
*dual -= target;
} else {
// Centering dual on primal > 0
*dual -= meanOnUnsaturated(*dual);
}
// Computing gradient norm
G = computeSquaredNorm(*dual);
// Updating search direction (conjugate gradient)
updateSearchDirection(delta * G / Gold);
Gold = G;
// Computing critical step
Real tau = computeCriticalStep(target);
// Update primal variable
delta = updatePrimal(tau);
// We should scale to constraint
if (constraint_type == variable_type)
*primal *= target / primal->mean();
error = computeError() / error_norm;
Real cost_f = functional.computeF(*primal, *dual);
printState(n, cost_f, error);
// printState(n, 0, error);
} while (error > this->tolerance && n++ < this->max_iterations);
// Final update of dual variable
functional.computeGradF(*primal, *dual);
// Enforcing dual constraint via gradient
if (constraint_type != variable_type) {
// Real dbar = meanOnUnsaturated(*dual);
// *dual += 2*target + dbar;
*dual -= target;
}
// Primal is pressure: need to make sure gap is positive
if (variable_type == pressure) {
*dual -= dual->min();
model.getDisplacement() = *dual;
model.getDisplacement() += this->surface;
}
// Primal is gap: need to make sure dual is positive (pressure + adhesion)
else {
model.getDisplacement() = *primal;
model.getDisplacement() += this->surface;
Real p_mean = -dual->min();
this->model.getBEEngine().solveDirichlet(model.getDisplacement(),
model.getTraction());
model.getTraction() += p_mean;
}
return error;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$
*/
Real PolonskyKeerRey::meanOnUnsaturated(const GridBase<Real>& field) const {
// Sum on unsaturated
Real sum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
*primal, field);
// Number of unsaturated points
UInt n_unsat = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal);
return sum / n_unsat;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
*/
Real PolonskyKeerRey::computeSquaredNorm(const GridBase<Real>& field) {
return Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; },
*primal, field);
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i'
* t_i} } \f$
*/
Real PolonskyKeerRey::computeCriticalStep(Real target) {
switch (variable_type) {
case gap:
model.getBEEngine().solveDirichlet(*search_direction,
*projected_search_direction);
break;
case pressure:
model.getBEEngine().solveNeumann(*search_direction,
*projected_search_direction);
break;
}
if (constraint_type == variable_type)
*projected_search_direction -=
meanOnUnsaturated(*projected_search_direction);
else {
// Real rbar = meanOnUnsaturated(*projected_search_direction);
// *projected_search_direction += 2 * target + rbar;
*projected_search_direction -= target;
}
Real num = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
return (p > 0) ? q * t : 0;
},
*primal, *dual, *search_direction);
Real denum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, Real & r, Real & t) {
return (p > 0) ? r * t : 0;
},
*primal, *projected_search_direction, *search_direction);
return num / denum;
}
/* -------------------------------------------------------------------------- */
/**
* Update steps:
* 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$
* 2. Truncate all \f$p\f$ negative
* 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i =
* p_i - \tau q_i\f$
*/
bool PolonskyKeerRey::updatePrimal(Real step) {
UInt na_num = Loop::reduce<operation::plus>(
[step] CUDA_LAMBDA(Real & p, Real & q, Real & t) -> UInt {
p -= step * t; // Updating primal
if (p < 0)
p = 0; // Truncating neg values
if (p == 0 && q < 0) { // If non-admissible state
p -= step * q;
return 1;
} else
return 0;
},
*primal, *dual, *search_direction);
return na_num == 0;
}
/* -------------------------------------------------------------------------- */
/**
* Error is based on \f$ \sum{p_i q_i} \f$
*/
Real PolonskyKeerRey::computeError() {
/// Making sure dual respects constraint
*dual -= dual->min();
Real error = primal->dot(*dual);
Real norm = 1;
if (variable_type == pressure)
norm = std::abs(primal->sum() * this->surface_stddev);
else
norm = std::abs(dual->sum() * this->surface_stddev);
norm *= primal->getNbPoints();
return std::abs(error) / norm;
}
/* -------------------------------------------------------------------------- */
/**
* Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t}
* \f$
*/
void PolonskyKeerRey::updateSearchDirection(Real factor) {
Loop::loop(
[factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
t = (p > 0) ? q + factor * t : 0;
},
*primal, *dual, *search_direction);
}
/* -------------------------------------------------------------------------- */
PolonskyKeerRey::~PolonskyKeerRey() {
switch (variable_type) {
case gap:
delete primal;
break;
case pressure:
delete dual;
break;
}
delete search_direction;
delete projected_search_direction;
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help