Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F59753520
polonsky_keer_rey.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Apr 24, 23:39
Size
9 KB
Mime Type
text/x-c
Expires
Fri, Apr 26, 23:39 (2 d)
Engine
blob
Format
Raw Data
Handle
17209130
Attached To
rTAMAAS tamaas
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 "fft_plan_manager.hh"
#include <iomanip>
#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
PolonskyKeerRey::PolonskyKeerRey(Model& model, const 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(new functional::ElasticFunctionalGap(
model.getBEEngine(), this->surface),
true);
break;
}
case pressure: {
primal = pressure_;
dual = gap_;
functional.addFunctionalTerm(new functional::ElasticFunctionalPressure(
model.getBEEngine(), this->surface),
true);
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 *= target / primal->mean();
else if (constraint_type != variable_type)
*primal = this->surface_stddev;
do {
// Computing functional gradient
functional.computeGradF(*primal, *dual);
Real dbar = meanOnUnsaturated(*dual);
// Enforcing dual constraint via gradient
if (constraint_type != variable_type) {
*dual += 2 * target + dbar;
} else {
// Centering dual on primal > 0
*dual -= dbar;
}
// 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);
} while (error > this->tolerance && n++ < this->max_iterations);
// Final update of dual variable
functional.computeGradF(*primal, *dual);
Real dual_min = dual->min();
*dual -= dual_min;
// Primal is pressure: need to make sure gap is positive
if (variable_type == pressure) {
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;
this->model.solveDirichlet();
model.getTraction() -= dual_min;
}
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) const {
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;
}
Real rbar = meanOnUnsaturated(*projected_search_direction);
if (constraint_type == variable_type)
*projected_search_direction -= rbar;
else {
*projected_search_direction += 2 * target + rbar;
}
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;
FFTPlanManager::get().clean(); // TODO find a smarter way to deal with plans
}
/* -------------------------------------------------------------------------- */
void PolonskyKeerRey::addFunctionalTerm(functional::Functional * func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment