Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/solvers/condat.cpp
b0c44a420f99archive/legacy_mpi

Last Change
Raw File

condat.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2024 Lucas Frérot
*
* 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 "condat.hh"
#include "logger.hh"
#include <iomanip>
/* -------------------------------------------------------------------------- */
namespace tamaas {
Condat::Condat(Model& model, const GridBase<Real>& surface, Real tolerance,
Real mu)
: Kato(model, surface, tolerance, mu) {
pressure_old =
allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
model.getTraction().getNbComponents());
}
/* -------------------------------------------------------------------------- */
Real Condat::solve(std::vector<Real> p0_vec) {
GridBase<Real> p0(p0_vec.size(), 1);
std::copy(p0_vec.begin(), p0_vec.end(), p0.begin());
TAMAAS_ASSERT(
p0.getNbPoints() == pressure->getNbComponents(),
"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, grad_step);
break;
case model_type::surface_2d:
cost = solveTmpl<model_type::surface_2d>(p0, grad_step);
break;
default:
break;
}
return cost;
}
template <model_type type>
Real Condat::solveTmpl(GridBase<Real>& p0, Real grad_step) {
// 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->uniformSetComponents(p0);
*pressure_old = 0;
Grid<Real, 1> q({1}, comp);
q = 0;
Real lipschitz = engine.getNeumannNorm();
Real sigma = 2 * grad_step / lipschitz;
Vector<Real, comp> h_mean = computeMean<comp>(*surfaceComp);
*surfaceComp -= h_mean;
do {
updateGap<comp>(sigma, grad_step, q);
*pressure -= *gap;
enforcePressureCoulomb<comp>();
updateLagrange<comp>(q, p0);
*pressure_old = *pressure;
cost = computeCost();
printState(n, cost, cost);
} while ((cost > this->tolerance || cost == 0.0) &&
n++ < this->max_iterations);
*surfaceComp += h_mean;
computeFinalGap<comp>();
return cost;
}
/* -------------------------------------------------------------------------- */
// template <UInt comp>
// void Condat::updateGap(Real sigma, Real grad_step, GridBase<Real>& q) {
// engine.solveNeumann(*pressure, *gap);
// VectorProxy<Real, comp> _q = q(0);
// Vector<Real, comp> qt = _q;
// qt *= (1 - grad_step);
// Loop::stridedLoop(
// [sigma, qt] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
// VectorProxy<Real, comp>&& h) {
// g -= h;
// g *= sigma;
// g += qt;
// },
// *gap, *surfaceComp);
// }
template <UInt comp>
void Condat::updateGap(Real sigma, Real grad_step, GridBase<Real>& q) {
computeGradient<comp>();
VectorProxy<Real, comp> _q(q(0));
Vector<Real, comp> qt = _q;
qt *= (1 - grad_step);
Loop::loop(
[sigma, qt] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
g *= sigma;
g += qt;
},
range<VectorProxy<Real, comp>>(*gap));
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Condat::updateLagrange(GridBase<Real>& q, GridBase<Real>& p0) {
*pressure_old *= -0.5;
*pressure_old += *pressure;
Vector<Real, comp> corr = computeMean<comp>(*pressure_old);
corr *= 2;
q += corr;
q -= p0;
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
c4science · Help