Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F60221866
accumulator.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Apr 28, 10:25
Size
7 KB
Mime Type
text/x-c++
Expires
Tue, Apr 30, 10:25 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17319838
Attached To
rTAMAAS tamaas
accumulator.hh
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/>.
*
*/
/* -------------------------------------------------------------------------- */
#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;
}
auto r = this->r, xc = this->xc; // for capture by value
Loop::loop(
[r, xc] CUDA_LAMBDA(VectorProxy<const Real, trait::dimension - 1> qv,
Source wk0, Source wk1, Source acc_g0,
Source 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
Event Timeline
Log In to Comment