Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/model/kelvin.cpp
7767124f567amaxwell_viscoelastic…

Last Change
Raw File

kelvin.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 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 "kelvin.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
Kelvin<type, derivative>::Kelvin(Model* model) : VolumePotential<type>(model) {
setIntegrationMethod(integration_method::linear, 1e-12);
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::setIntegrationMethod(integration_method method,
Real cutoff) {
this->method = method;
this->cutoff = cutoff;
Logger logger;
if (this->method == integration_method::linear) {
logger.get(LogLevel::debug)
<< TAMAAS_MSG("Setting linear integration method");
this->initialize(dtrait::template source_components<type>,
dtrait::template out_components<type>,
this->model->getDiscretization()[0]);
}
else {
logger.get(LogLevel::debug) << TAMAAS_MSG(
"Setting cutoff integration method (cutoff ", this->cutoff, ')');
this->initialize(dtrait::template source_components<type>,
dtrait::template out_components<type>, 1);
}
auto max_q = Loop::reduce<operation::max>(
[] CUDA_LAMBDA(VectorProxy<const Real, trait::boundary_dimension> qv) {
return qv.l2norm();
},
range<VectorProxy<const Real, trait::boundary_dimension>>(
this->wavevectors));
if (this->method == integration_method::linear and
not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
logger.get(LogLevel::warning)
<< "Probable overflow of integral computation (consider "
"changing integration method to integration_method::cutoff or "
"compiling with real_type='long double')\n";
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::applyIf(GridBase<Real>& source,
GridBase<Real>& out,
filter_t pred) const {
KelvinInfluence kelvin(this->model->getShearModulus(),
this->model->getPoissonRatio());
parent::transformSource(source, pred);
// Reset buffer values
for (auto&& layer : this->out_buffer)
layer = 0;
if (method == integration_method::linear) {
linearIntegral(out, kelvin);
} else {
cutoffIntegral(out, kelvin);
}
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::linearIntegral(GridBase<Real>& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper<type, KelvinInfluence> helper;
helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
this->model->getSystemSize().front(), kelvin);
helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
helper.makeFundamentalGreatAgain(this->out_buffer);
parent::transformOutput(
[](auto&& out_buffer, auto layer) ->
typename parent::BufferType& { return out_buffer[layer]; },
out);
}
template <model_type type, UInt derivative>
void Kelvin<type, derivative>::cutoffIntegral(GridBase<Real>& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper<type, KelvinInfluence> helper;
auto func = [&](auto&& out_buffer,
auto layer) -> typename parent::BufferType& {
auto&& out = out_buffer.front();
out = 0;
helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors,
this->model->getSystemSize().front(), cutoff, kelvin);
helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
helper.makeFundamentalGreatAgain(out);
return out;
};
parent::transformOutput(func, out);
}
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
std::pair<UInt, UInt> Kelvin<type, derivative>::matvecShape() const {
const auto N = this->model->getDisplacement().getGlobalNbPoints();
return std::make_pair(dtrait::template out_components<type> * N,
dtrait::template source_components<type> * N);
}
/* -------------------------------------------------------------------------- */
template <model_type type, UInt derivative>
GridBase<Real> Kelvin<type, derivative>::matvec(GridBase<Real>& X) const {
constexpr auto dim = trait::dimension;
const auto shape = this->model->getDiscretization();
Grid<Real, dim> x{shape, dtrait::template source_components<type>, X.view()};
Grid<Real, dim> y{shape, dtrait::template out_components<type>};
this->apply(x, y);
return y;
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Kelvin<model_type::volume_2d, 0>;
template class Kelvin<model_type::volume_2d, 1>;
template class Kelvin<model_type::volume_2d, 2>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
c4science · Help