Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/model/integration/accumulator.hh
c74466647365features/declarative

Last Change
Raw File

accumulator.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef ACCUMULATOR_HH
#define ACCUMULATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid_hermitian.hh"
#include "integrator.hh"
#include "model_type.hh"
#include "static_types.hh"
#include <array>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Accumulation integral manager
template <model_type type, typename Source>
class Accumulator {
using trait = model_type_traits<type>;
using BufferType = GridHermitian<Real, trait::boundary_dimension>;
public:
/// Constructor
Accumulator(
std::vector<BufferType>& node_values,
const std::vector<UInt>& hermitian_sizes,
const Grid<Real, trait::boundary_dimension>& wavevectors)
: node_values(node_values), wavevectors(wavevectors) {
// Resizing accumulators (with source)
std::for_each(this->accumulator.begin(), this->accumulator.end(),
[&](auto& accumulator) {
accumulator.setNbComponents(Source::size);
accumulator.resize(hermitian_sizes);
});
}
/// Initialize uniform mesh
void makeUniformMesh(Real domain_size) {
UInt N = node_values.size();
Real dx = domain_size / (N - 1);
Real x = 0;
node_positions.resize(N);
for (UInt i = 0; i < N; ++i, x += dx) {
node_positions[i] = x;
}
}
/// Direction flag
enum class direction { forward, backward };
/// direction <=> upper
static constexpr bool is_upper(direction dir) {
return dir == direction::backward;
}
/// Forward/Backward iterator for integration passes
template <direction dir>
class iterator {
public:
/// Constructor
iterator(Accumulator& acc, Int l) : acc(acc), l(l) {}
/// Copy constructor
iterator(const iterator& o) : acc(o.acc), l(o.l) {}
/// Compare
bool operator!=(const iterator& o) const { return l != o.l; }
/// Increment
iterator& operator++() {
// current layer index
l += (dir == direction::forward) ? 1 : -1;
// part of the domain to be integrated
constexpr bool upper = is_upper(dir);
// value of k to be added to the sum
const auto k = (dir == direction::forward) ? l - 1 : l;
const auto r = 0.5 * (acc.node_positions[k + 1] - acc.node_positions[k]);
const auto xc = 0.5 * (acc.node_positions[k + 1] + acc.node_positions[k]);
using integ = Integrator<1>;
Loop::loop(
[&](auto qv, auto wk0, auto wk1, auto acc_g0, auto acc_g1) {
Real q = qv.l2norm();
acc_g0 += wk0 * integ::template G0<upper, 0>(q, r, xc) +
wk1 * integ::template G0<upper, 1>(q, r, xc);
acc_g1 += wk0 * integ::template G1<upper, 0>(q, r, xc) +
wk1 * integ::template G1<upper, 1>(q, r, xc);
},
range<VectorProxy<const Real, trait::dimension - 1>>(acc.wavevectors),
range<Source>(acc.node_values[k]),
range<Source>(acc.node_values[k + 1]),
range<Source>(acc.accumulator.front()),
range<Source>(acc.accumulator.back()));
return *this;
}
/// Dereference iterator
std::tuple<int, Real, BufferType&, BufferType&> operator*() {
return {l, acc.node_positions[l], acc.accumulator.front(),
acc.accumulator.back()};
}
private:
Accumulator& acc;
Int l;
};
/// Range for convenience
template <direction dir>
struct acc_range {
iterator<dir> _begin, _end;
iterator<dir> begin() const { return _begin; }
iterator<dir> end() const { return _end; }
};
/// Prepare forward loop
acc_range<direction::forward> forward() {
reset();
return acc_range<direction::forward>{
iterator<direction::forward>(*this, 0),
iterator<direction::forward>(*this, node_positions.size())};
}
/// Prepare backward loop
acc_range<direction::backward> backward() {
reset();
return acc_range<direction::backward>{
iterator<direction::backward>(*this, node_positions.size() - 1),
iterator<direction::backward>(*this, -1)};
}
void reset() {
for (auto&& acc : accumulator)
acc = 0;
}
private:
std::array<BufferType, 2> accumulator;
std::vector<Real> node_positions;
std::vector<BufferType>& node_values;
const Grid<Real, trait::boundary_dimension>& wavevectors;
};
} // namespace tamaas
#endif // ACCUMULATOR_HH
c4science · Help