Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/condat.cpp
075c8d026736features/netcdf-mpi

Last Change
Raw File

condat.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 "condat.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Condat::Condat(Model& model, const GridBase<Real>& surface,
Real tolerance, Real mu)
: Kato(model, surface, tolerance, mu) {}
/* -------------------------------------------------------------------------- */
Real Condat::solve(GridBase<Real>& p0) {
if (p0.getNbPoints() != comp) {
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 Condat::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;
Real cost = 0;
UInt n = 0;
pressure->uniformSetComponents(p0);
GridBase<Real> pressure_old(*pressure); pressure_old = 0;
GridBase<Real> q(*pressure); q = 0;
Real lipschitz = engine.getNeumannNorm();
Real gamma = 0.9;
Real sigma = 2 * gamma / lipschitz;
do {
updateGap<type>(sigma, gamma, q);
*pressure -= *gap;
enforcePressureCoulomb<type>();
updateLagrange<type>(q, pressure_old, p0);
pressure_old = *pressure;
cost = computeCost();
printState(n, cost, cost);
} while (cost > this->tolerance && n++ < this->max_iterations);
engine.solveNeumann(*pressure, *gap);
*gap -= surfaceComp;
Real beta = computeBeta<type>();
computeFinalGap<type>(beta);
model.getDisplacement() = *gap;
return cost;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Condat::updateGap(Real sigma, Real gamma, GridBase<Real>& q) {
constexpr UInt comp = model_type_traits<type>::components;
engine.solveNeumann(*pressure, *gap);
Loop::stridedLoop(
[sigma, gamma] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
VectorProxy<Real, comp>&& h,
VectorProxy<Real, comp>&& q) {
g -= h;
g *= sigma;
Vector<Real, comp> qt = q;
qt *= (1 - gamma);
g += qt;
},
*gap, surfaceComp, q);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void Condat::updateLagrange(GridBase<Real>& q, GridBase<Real>& p_old,
GridBase<Real>& p0) {
constexpr UInt comp = model_type_traits<type>::components;
// GridBase<Real> corr(p0);
Grid<Real, 1> corr({1}, comp);
p_old *= -0.5;
p_old += *pressure;
computeMean<comp>(p_old, corr);
corr *= 2;
corr -= p0;
addUniform<comp>(q, corr);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
c4science · Help