Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/model/integration/accumulator.hh
d5841e5c2852maxwell_viscoelastic…

Last Change
Raw File

accumulator.hh
View Options

/**
* @file
* LICENSE
*
* Copyright (©) 2016-2021 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef ACCUMULATOR_HH
#define ACCUMULATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid_hermitian.hh"
#include "integrator.hh"
#include "model_type.hh"
#include "static_types.hh"
#include <array>
#include <thrust/iterator/zip_iterator.h>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Accumulation integral manager
template <model_type type, typename Source,
typename = std::enable_if_t<is_proxy<Source>::value>>
class Accumulator {
using trait = model_type_traits<type>;
using BufferType = GridHermitian<Real, trait::boundary_dimension>;
public:
/// Constructor
Accumulator() {
// Resizing accumulators (with source)
for (auto&& acc : accumulator)
acc.setNbComponents(Source::size);
}
/// Initialize uniform mesh
void makeUniformMesh(UInt N, Real domain_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;
}
}
const std::vector<Real>& nodePositions() const { return node_positions; }
/// Direction flag
enum class direction { forward, backward };
/// Forward/Backward iterator for integration passes
template <direction dir>
struct iterator {
using integ = Integrator<1>;
/// Constructor
iterator(Accumulator& acc, Int k) : acc(acc), k(k) {
l = (dir != direction::forward) ? k + 1 : k;
}
/// Copy constructor
iterator(const iterator& o)
: acc(o.acc), k(o.k), l(o.l), r(o.r), xc(o.xc) {}
static constexpr bool upper = dir == direction::backward;
/// Compare
bool operator!=(const iterator& o) const { return k != o.k; }
/// Increment
iterator& operator++() {
// update current layer and element. skip last operation
if (not setup()) {
next();
return *this;
}
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);
acc_g0 += wk1 * integ::template G0<upper, 1>(q, r, xc);
acc_g1 += wk0 * integ::template G1<upper, 0>(q, r, xc);
acc_g1 += wk1 * integ::template G1<upper, 1>(q, r, xc);
},
range<VectorProxy<const Real, trait::dimension - 1>>(
acc.waveVectors())
.headless(),
range<Source>(acc.nodeValues()[k]).headless(),
range<Source>(acc.nodeValues()[k + 1]).headless(),
range<Source>(acc.accumulator.front()).headless(),
range<Source>(acc.accumulator.back()).headless());
if (mpi::rank() == 0) {
// accumulating the fundamental mode
Source wk0(&this->acc.nodeValues()[this->k](0)),
wk1(&this->acc.nodeValues()[this->k + 1](0));
Source acc_tensor(&acc.accumulator.front()(0));
acc_tensor += wk0 * integ::template F<0>(this->r);
acc_tensor += wk1 * integ::template F<1>(this->r);
}
next();
return *this;
}
/// Dereference iterator (TODO uniformize tuple types)
std::tuple<Int, Real, BufferType&, BufferType&> operator*() {
return {l, acc.node_positions[l], acc.accumulator.front(),
acc.accumulator.back()};
}
protected:
/// Update index layer and element info
bool setup() {
if (k + 1 >= static_cast<Int>(acc.node_positions.size()) or k < 0)
return false;
r = 0.5 * (acc.node_positions[k + 1] - acc.node_positions[k]);
xc = 0.5 * (acc.node_positions[k + 1] + acc.node_positions[k]);
return true;
}
/// Set current layer and update element index
void next() {
l = layer(k);
k = nextElement(k);
}
/// Element index => layer index
static Int layer(Int element) {
return element + (dir == direction::forward);
}
/// Next element index in right direction
static Int nextElement(Int element) {
return element + ((dir == direction::forward) ? 1 : -1);
}
protected:
Accumulator& acc; /// accumulator holder
Int k = 0; ///< element index
Int l = 0; ///< layer index
Real r = 0; ///< element radius
Real xc = 0; ///< element center
};
/// Range for convenience
template <typename Iterator>
struct acc_range {
Iterator _begin, _end;
Iterator begin() const { return _begin; }
Iterator end() const { return _end; }
};
/// Prepare forward loop
acc_range<iterator<direction::forward>>
forward(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
using it = iterator<direction::forward>;
reset(nvalues, wvectors);
return acc_range<it>{it(*this, 0), it(*this, node_positions.size())};
}
/// Prepare backward loop
acc_range<iterator<direction::backward>>
backward(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
using it = iterator<direction::backward>;
reset(nvalues, wvectors);
// don't ask why the range has those values
// if it works, don't change it
return acc_range<it>{it(*this, node_positions.size() - 2), it(*this, -2)};
}
void reset(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
node_values = &nvalues;
wavevectors = &wvectors;
for (auto&& acc : accumulator) {
acc.resize(wvectors.sizes());
acc = 0;
}
}
std::vector<BufferType>& nodeValues() {
TAMAAS_ASSERT(node_values, "Node values is invalid");
return *node_values;
}
const Grid<Real, trait::boundary_dimension>& waveVectors() {
TAMAAS_ASSERT(wavevectors, "Wavevectors is invalid");
return *wavevectors;
}
/// Create a range over the elements in the mesh
auto elements() {
auto range = Loop::range(node_positions.size() - 1);
auto begin = thrust::make_zip_iterator(thrust::make_tuple(
range.begin(), node_positions.begin(), node_positions.begin() + 1));
auto end = thrust::make_zip_iterator(thrust::make_tuple(
range.end(), node_positions.end() - 1, node_positions.end()));
return acc_range<decltype(begin)>{begin, end};
}
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