Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/solvers/kato_saturated.cpp
43f7d9ba9056maxwell_viscoelastic…

Last Change
Raw File

kato_saturated.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 "kato_saturated.hh"
#include "logger.hh"
#include <iomanip>
#include <limits>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
KatoSaturated::KatoSaturated(Model& model, const GridBase<Real>& surface,
Real tolerance, Real pmax)
: PolonskyKeerRey(model, surface, tolerance, PolonskyKeerRey::pressure,
PolonskyKeerRey::pressure),
pmax(pmax) {}
/* -------------------------------------------------------------------------- */
Real KatoSaturated::solve(std::vector<Real> load) {
GridBase<Real> initial_surface = surface;
GridBase<Real> residual_disp = surface;
residual_disp = 0;
auto negative = [this](const GridBase<Real>& f) {
return f.min() < -std::abs(this->tolerance);
};
UInt n = 0;
do {
surface = initial_surface;
surface += residual_disp;
PolonskyKeerRey::solve(load);
// Update the rough surface
Loop::loop([](Real& h_pl, Real g) { h_pl += g * (g < 0); }, residual_disp,
*this->dual);
} while (negative(*this->dual) and n++ < this->max_iterations);
surface = initial_surface;
*this->displacement_view -= residual_disp;
return Real(n >= this->max_iterations);
}
/* -------------------------------------------------------------------------- */
Real KatoSaturated::meanOnUnsaturated(const GridBase<Real>& /*field*/) const {
return 0;
}
Real KatoSaturated::computeSquaredNorm(const GridBase<Real>& /*field*/) const {
return 1.;
}
void KatoSaturated::updateSearchDirection(Real /*factor*/) {
*this->search_direction = *this->dual;
}
Real KatoSaturated::computeCriticalStep(Real /*target*/) {
// integral_op->apply(*search_direction, *projected_search_direction);
// Real num = search_direction->dot(*search_direction);
// Real denum = projected_search_direction->dot(*search_direction);
// return 0.1 * num / denum;
return 1;
}
bool KatoSaturated::updatePrimal(Real step) {
UInt na_num = Loop::reduce<operation::plus>(
[step] CUDA_LAMBDA(Real & p, const Real& /*q*/, const Real& t) -> UInt {
p -= step * t; // Updating primal
return 0;
// 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;
}
/* -------------------------------------------------------------------------- */
Real KatoSaturated::computeError() {
// We shift the gap by the minimum on unsaturated area
Real shift = Loop::reduce<operation::min>(
[&] CUDA_LAMBDA(const auto& p, const auto& d) {
return (p < pmax) ? d : std::numeric_limits<Real>::max();
},
*primal, *dual);
*dual -= shift;
Real norm = 1;
// Ignore points that are saturated
Real error = Loop::reduce<operation::plus>(
[&] CUDA_LAMBDA(auto p, auto d) { return (p < pmax) ? p * d : 0; },
*primal, *dual);
if (std::isnan(error))
TAMAAS_EXCEPTION("Encountered NaN in complementarity error: this may be "
"caused by a contact area of a single node.");
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;
}
/* -------------------------------------------------------------------------- */
void KatoSaturated::enforceMeanValue(Real mean) {
// We want to cancel the difference between saturated alpha + t and the
// applied pressure
*primal -= primal->mean();
auto f = [&](Real scale) {
Real sum = Loop::reduce<operation::plus>(
[&] CUDA_LAMBDA(Real t) -> Real {
t += scale;
if (t > pmax)
return pmax;
else if (t < 0)
return 0;
else
return t;
},
*this->primal);
sum /= this->primal->getNbPoints();
return sum - mean;
};
if (pmax < mean)
TAMAAS_EXCEPTION("Cannot find equilibrium");
// Dichotomy + Secant Newton on f
// Initial points
Real x_n_2 = -primal->max(), x_n_1 = -primal->min(), x_n = 0.;
Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
// Find a not-to-large search interval
Logger().get(LogLevel::debug) << "Searching equilibrium interval\n";
while (std::signbit(f(x_n_1)) == std::signbit(f(x_n_2))) {
x_n_1 *= 10;
x_n_2 *= 10;
}
Logger().get(LogLevel::debug)
<< "Reducing interval [abs(x1/x2) = " << std::abs(x_n_1 / x_n_2) << "]\n";
UInt n_dic = 10;
f_n_1 = f(x_n_1);
for (UInt i = 0; i < n_dic; ++i) {
x_n = 0.5 * (x_n_1 + x_n_2);
f_n = f(x_n);
if (std::signbit(f_n) == std::signbit(f_n_1)) {
x_n_1 = x_n;
f_n_1 = f_n;
} else
x_n_2 = x_n;
}
Logger().get(LogLevel::debug)
<< "Starting Newton secant [abs(x1/x2) = " << std::abs(x_n_1 / x_n_2)
<< "]\n";
// Secant loop
do {
f_n_2 = f(x_n_2);
f_n_1 = f(x_n_1);
if (f_n_1 == f_n_2) break; // Avoid nans
x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
f_n = f(x_n);
x_n_2 = x_n_1;
x_n_1 = x_n;
} while (std::abs(f_n / mean) > 1e-14);
// Pressure update
Loop::loop(
[&] CUDA_LAMBDA(auto& t) {
t += x_n;
if (t > pmax)
t = pmax;
else if (t < 0)
t = 0.;
},
*this->primal);
}
/* -------------------------------------------------------------------------- */
void KatoSaturated::enforceAdmissibleState() {
/// Make dual admissible
Real shift = Loop::reduce<operation::min>(
[&] CUDA_LAMBDA(const auto& p, const auto& d) {
return (p < pmax) ? d : std::numeric_limits<Real>::max();
},
*primal, *dual);
*dual -= shift;
*displacement_view = *dual;
*displacement_view += this->surface;
}
} // namespace tamaas
c4science · Help