Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/beck_teboulle.cpp
a240eb50bf7efeatures/netcdf-mpi

Last Change
Raw File

beck_teboulle.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 "beck_teboulle.hh"
#include "logger.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BeckTeboulle::BeckTeboulle(Model& model, const GridBase<Real>& surface,
Real tolerance, Real mu)
: Kato(model, surface, tolerance, mu) {}
/* -------------------------------------------------------------------------- */
Real BeckTeboulle::solve(GridBase<Real>& g0) {
if (g0.getNbPoints() != pressure->getNbComponents()) {
TAMAAS_EXCEPTION(
"Target mean gap 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>(g0);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(g0);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real BeckTeboulle::solveTmpl(GridBase<Real>& g0) {
// Printing column headers
Logger().get(LogLevel::info) << 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;
UInt n = 0;
*pressure = 0;
GridBase<Real> z(*pressure);
GridBase<Real> pressure_old(*pressure);
Real t0 = 1;
Real t1 = 1;
Real m = 0;
Real lipschitz = engine.getNeumannNorm();
do {
engine.solveNeumann(*pressure, *gap);
addUniform<comp>(*gap, g0);
*gap -= *surfaceComp;
z = *gap;
z /= -lipschitz;
*pressure += z;
enforcePressureCoulomb<comp>();
z = *pressure;
t1 = 1 + std::sqrt(4 * t0 * t0 + 1) / 2;
m = 1 + (t0 - 1) / t1;
z -= pressure_old;
z *= m;
z += pressure_old;
t0 = t1;
pressure_old = *pressure;
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
computeFinalGap<comp>();
return cost;
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help