Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/core/grid_base.hh
1fc7644bad52maxwell_viscoelastic…

Last Change
Raw File

grid_base.hh
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* 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 GRID_BASE_HH
#define GRID_BASE_HH
/* -------------------------------------------------------------------------- */
#include "array.hh"
#include "iterator.hh"
#include "loop.hh"
#include "mpi_interface.hh"
#include "static_types.hh"
#include "tamaas.hh"
#include <cstddef>
#include <limits>
#include <utility>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Dimension agnostic grid with components stored per points
*/
template <typename T>
class GridBase {
public:
/// Constructor by default
GridBase() = default;
/// Copy constructor
GridBase(const GridBase& o) { this->copy(o); }
/// Move constructor (transfers data ownership)
GridBase(GridBase&& o) noexcept
: data(std::move(o.data)),
nb_components(std::exchange(o.nb_components, 1)) {}
/// Initialize with size
explicit GridBase(UInt nb_points, UInt nb_components = 1) {
setNbComponents(nb_components);
this->resize(nb_points * nb_components);
}
/// Destructor
virtual ~GridBase() = default;
/* ------------------------------------------------------------------------ */
/* Iterator class */
/* ------------------------------------------------------------------------ */
public:
using iterator = iterator_::iterator<T>;
using const_iterator = iterator_::iterator<const T>;
using value_type = T;
using reference = value_type&;
public:
/// Get grid dimension
virtual UInt getDimension() const { return 1; }
/// Get internal data pointer (const)
const T* getInternalData() const { return this->data.data(); }
/// Get internal data pointer (non-const)
T* getInternalData() { return this->data.data(); }
/// Get number of components
UInt getNbComponents() const { return nb_components; }
/// Set number of components
void setNbComponents(UInt n) { nb_components = n; }
/// Get total size
virtual UInt dataSize() const { return this->data.size(); }
/// Get global data size
UInt globalDataSize() const {
return mpi::allreduce<operation::plus>(dataSize());
}
/// Get number of points
UInt getNbPoints() const {
return this->dataSize() / this->getNbComponents();
}
/// Get global number of points
UInt getGlobalNbPoints() const {
return this->globalDataSize() / this->getNbComponents();
}
/// Set components
void uniformSetComponents(const GridBase<T>& vec);
/// Resize
void resize(UInt size) { this->data.assign(size, T(0.)); }
/// Reserve storage w/o changing data logic
void reserve(UInt size) { this->data.reserve(size); }
/* ------------------------------------------------------------------------ */
/* Iterator methods */
/* ------------------------------------------------------------------------ */
virtual iterator begin(UInt n = 1) {
return iterator(this->getInternalData(), n);
}
virtual iterator end(UInt n = 1) {
return iterator(this->getInternalData() + this->dataSize(), n);
}
virtual const_iterator begin(UInt n = 1) const {
return const_iterator(this->getInternalData(), n);
}
virtual const_iterator end(UInt n = 1) const {
return const_iterator(this->getInternalData() + this->dataSize(), n);
}
/* ------------------------------------------------------------------------ */
/* Operators */
/* ------------------------------------------------------------------------ */
inline T& operator()(UInt i) { return this->data[i]; }
inline const T& operator()(UInt i) const { return this->data[i]; }
#define VEC_OPERATOR(op) \
template <typename T1> \
GridBase& operator op(const GridBase<T1>& other)
VEC_OPERATOR(+=);
VEC_OPERATOR(*=);
VEC_OPERATOR(-=);
VEC_OPERATOR(/=);
#undef VEC_OPERATOR
template <typename T1>
T dot(const GridBase<T1>& other) const;
#define SCALAR_OPERATOR(op) GridBase& operator op(const T& e)
SCALAR_OPERATOR(+=);
SCALAR_OPERATOR(*=);
SCALAR_OPERATOR(-=);
SCALAR_OPERATOR(/=);
SCALAR_OPERATOR(=);
#undef SCALAR_OPERATOR
/// Broadcast operators
#define BROADCAST_OPERATOR(op) \
template <typename DT, typename ST, UInt size> \
GridBase& operator op(const StaticArray<DT, ST, size>& b)
BROADCAST_OPERATOR(+=);
BROADCAST_OPERATOR(-=);
BROADCAST_OPERATOR(*=);
BROADCAST_OPERATOR(/=);
BROADCAST_OPERATOR(=);
#undef BROADCAST_OPERATOR
/* ------------------------------------------------------------------------ */
/* Min/max */
/* ------------------------------------------------------------------------ */
T min() const;
T max() const;
T sum() const;
T mean() const { return sum() / static_cast<T>(globalDataSize()); }
T var() const;
T l2norm() const { return std::sqrt(dot(*this)); }
/* ------------------------------------------------------------------------ */
/* Move/Copy */
/* ------------------------------------------------------------------------ */
public:
GridBase& operator=(const GridBase& o) {
this->copy(o);
return *this;
}
GridBase& operator=(GridBase& o) {
this->copy(o);
return *this;
}
GridBase& operator=(GridBase&& o) noexcept {
this->move(std::move(o));
return *this;
}
template <typename T1>
void copy(const GridBase<T1>& other) {
if (other.dataSize() != this->dataSize())
this->resize(other.dataSize());
thrust::copy(other.begin(), other.end(), this->begin());
nb_components = other.nb_components;
}
template <typename T1>
void move(GridBase<T1>&& other) noexcept {
data = std::move(other.data);
nb_components = std::exchange(other.nb_components, 1);
}
GridBase& wrap(GridBase& o) {
data.wrap(o.data);
nb_components = o.nb_components;
return *this;
}
GridBase& wrap(const GridBase& o) {
data.wrap(o.data);
nb_components = o.nb_components;
return *this;
}
protected:
Array<T> data;
UInt nb_components = 1;
};
/* -------------------------------------------------------------------------- */
template <typename T>
void GridBase<T>::uniformSetComponents(const GridBase<T>& vec) {
if (vec.dataSize() != this->nb_components)
TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
auto begin_it{begin()};
auto end_it{end()}; // silencing nvcc warnings
for (auto it = begin(); it < end_it; ++it) {
UInt i = it - begin_it;
*it = vec(i % this->nb_components);
}
}
/* --------------------------------------------------------------------------
*/
#define SCALAR_OPERATOR_IMPL(op) \
template <typename T> \
inline GridBase<T>& GridBase<T>::operator op(const T& e) { \
Loop::loop([e] CUDA_LAMBDA(T& val) { val op e; }, *this); \
return *this; \
}
SCALAR_OPERATOR_IMPL(+=)
SCALAR_OPERATOR_IMPL(*=)
SCALAR_OPERATOR_IMPL(-=)
SCALAR_OPERATOR_IMPL(/=)
SCALAR_OPERATOR_IMPL(=)
#undef SCALAR_OPERATOR_IMPL
/* -------------------------------------------------------------------------- */
#define VEC_OPERATOR_IMPL(op) \
template <typename T> \
template <typename T1> \
inline GridBase<T>& GridBase<T>::operator op(const GridBase<T1>& other) { \
TAMAAS_ASSERT(other.dataSize() == this->dataSize(), \
"surface size does not match"); \
Loop::loop([] CUDA_LAMBDA(T& x, const T1& y) { x op y; }, *this, other); \
return *this; \
}
VEC_OPERATOR_IMPL(+=)
VEC_OPERATOR_IMPL(-=)
VEC_OPERATOR_IMPL(*=)
VEC_OPERATOR_IMPL(/=)
#undef VEC_OPERATOR
/* -------------------------------------------------------------------------- */
#define BROADCAST_OPERATOR_IMPL(op) \
template <typename T> \
template <typename DT, typename ST, UInt size> \
GridBase<T>& GridBase<T>::operator op(const StaticArray<DT, ST, size>& b) { \
TAMAAS_ASSERT(this->dataSize() % size == 0, \
"Broadcast operator cannot broadcast" \
<< size << " on array of size " << this->dataSize()); \
Tensor<StaticArray, typename std::remove_const<DT>::type, size> a(b); \
Loop::loop([a] CUDA_LAMBDA(VectorProxy<T, size> g) { g op a; }, \
range<VectorProxy<T, size>>(*this)); \
return *this; \
}
BROADCAST_OPERATOR_IMPL(+=)
BROADCAST_OPERATOR_IMPL(-=)
BROADCAST_OPERATOR_IMPL(*=)
BROADCAST_OPERATOR_IMPL(/=)
BROADCAST_OPERATOR_IMPL(=)
#undef BROADCAST_OPERATOR_IMPL
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::min() const {
// const auto id = [] CUDA_LAMBDA (const T&x){return x;};
return Loop::reduce<operation::min>([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::max() const {
return Loop::reduce<operation::max>([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::sum() const {
return Loop::reduce<operation::plus>([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T GridBase<T>::var() const {
const T mu = mean();
const T var = Loop::reduce<operation::plus>(
[mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
return var / (globalDataSize() - 1);
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename T1>
T GridBase<T>::dot(const GridBase<T1>& other) const {
return Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
}
} // namespace tamaas
#endif // GRID_BASE_HH
c4science · Help