diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 649ec51..4ac0e20 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,202 +1,202 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include "regularized_powerlaw.hh"
#include "surface_generator_filter.hh"
#include "surface_generator_filter_fft.hh"
#include "surface_generator_random_phase.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
class PyFilter : public Filter {
public:
using Filter::Filter;
// Overriding pure virtual functions
void
computeFilter(GridHermitian& filter_coefficients) const override {
PYBIND11_OVERLOAD_PURE(void, Filter, computeFilter,
filter_coefficients);
}
};
template
void wrapFilter(py::module& mod) {
auto name = makeDimensionName("Filter", dim);
py::class_, std::shared_ptr>, PyFilter>(
mod, name.c_str())
.def(py::init<>())
.def("computeFilter",
(void (Filter::*)(GridHermitian&) const) &
Filter::computeFilter);
}
/* -------------------------------------------------------------------------- */
template
void wrapIsopowerlaw(py::module& mod) {
std::string name = makeDimensionName("Isopowerlaw", dim);
py::class_, Filter, std::shared_ptr>>(
mod, name.c_str())
.def(py::init<>())
.def_property("q0", &Isopowerlaw::getQ0, &Isopowerlaw::setQ0)
.def_property("q1", &Isopowerlaw::getQ1, &Isopowerlaw::setQ1)
.def_property("q2", &Isopowerlaw::getQ2, &Isopowerlaw::setQ2)
.def_property("hurst", &Isopowerlaw::getHurst,
&Isopowerlaw::setHurst)
.def("rmsHeights", &Isopowerlaw::rmsHeights)
.def("moments", &Isopowerlaw::moments)
.def("alpha", &Isopowerlaw::alpha)
.def("rmsSlopes", &Isopowerlaw::rmsSlopes)
// legacy wrapper code
.def("getQ0", &Isopowerlaw::getQ0,
py::return_value_policy::reference)
.def("getQ1", &Isopowerlaw::getQ1,
py::return_value_policy::reference)
.def("getQ2", &Isopowerlaw::getQ2,
py::return_value_policy::reference)
.def("setQ0", &Isopowerlaw::setQ0,
py::return_value_policy::reference)
.def("setQ1", &Isopowerlaw::setQ1,
py::return_value_policy::reference)
.def("setQ2", &Isopowerlaw::setQ2,
py::return_value_policy::reference)
.def("setHurst", &Isopowerlaw::setHurst)
.def("getHurst", &Isopowerlaw::getHurst,
py::return_value_policy::reference);
name = makeDimensionName("RegularizedPowerlaw", dim);
py::class_, Filter,
std::shared_ptr>>(mod, name.c_str())
.def(py::init<>())
.def_property("q1", &RegularizedPowerlaw::getQ1,
&RegularizedPowerlaw::setQ1)
.def_property("q2", &RegularizedPowerlaw::getQ2,
&RegularizedPowerlaw::setQ2)
.def_property("hurst", &RegularizedPowerlaw::getHurst,
&RegularizedPowerlaw::setHurst);
}
/* -------------------------------------------------------------------------- */
template
void wrapSurfaceGenerators(py::module& mod) {
std::string generator_name = makeDimensionName("SurfaceGenerator", dim);
py::class_>(mod, generator_name.c_str())
.def("buildSurface", &SurfaceGenerator::buildSurface,
py::return_value_policy::reference_internal)
.def("setSizes",
- py::overload_cast&>(
+ py::overload_cast>(
&SurfaceGenerator::setSizes),
"surface_sizes"_a)
.def_property("random_seed", &SurfaceGenerator::getRandomSeed,
&SurfaceGenerator::setRandomSeed);
std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim);
py::class_, SurfaceGenerator>(
mod, filter_name.c_str())
.def(py::init<>())
.def("setFilter", &SurfaceGeneratorFilter::setFilter,
"Set PSD filter", "filter"_a)
.def("setSpectrum", &SurfaceGeneratorFilter::setSpectrum,
"Set PSD spectrum", "filter"_a)
// legacy wrapper code
.def("setRandomSeed", &SurfaceGeneratorFilter::setRandomSeed,
"[[Deprecated: use property]]");
std::string random_name =
makeDimensionName("SurfaceGeneratorRandomPhase", dim);
py::class_, SurfaceGenerator>(
mod, random_name.c_str())
.def(py::init<>())
.def("setFilter", &SurfaceGeneratorRandomPhase::setFilter,
"Set PSD filter", "filter"_a)
.def("setSpectrum", &SurfaceGeneratorRandomPhase::setSpectrum,
"Set PSD Spectrum", "spectrum"_a);
}
/* -------------------------------------------------------------------------- */
// Legacy wrap
void wrapSurfaceGeneratorFilterFFT(py::module& mod) {
py::class_(mod, "SurfaceGeneratorFilterFFT")
.def(py::init<>())
.def("getQ0",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getQ0());
})
.def("getQ1",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getQ1());
})
.def("getQ2",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getQ2());
})
.def("getHurst",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getHurst());
})
.def("analyticalRMS", &SurfaceGeneratorFilterFFT::analyticalRMS)
.def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface,
py::return_value_policy::reference_internal)
.def("getGridSize",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getGridSize());
})
.def("getRandomSeed",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getRandomSeed());
})
.def("getRMS",
[](SurfaceGeneratorFilterFFT& f) {
return smart_pointer(&f.getRMS());
})
.def("Init", &SurfaceGeneratorFilterFFT::Init);
}
/* -------------------------------------------------------------------------- */
void wrapSurface(py::module& mod) {
wrapFilter<1>(mod);
wrapFilter<2>(mod);
wrapIsopowerlaw<1>(mod);
wrapIsopowerlaw<2>(mod);
wrapSurfaceGenerators<1>(mod);
wrapSurfaceGenerators<2>(mod);
// legacy wrap
#if defined(LEGACY_BEM)
wrapSurfaceGeneratorFilterFFT(mod);
#endif
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/core/fftw_interface.hh b/src/core/fftw_interface.hh
index d2f9767..3be0603 100644
--- a/src/core/fftw_interface.hh
+++ b/src/core/fftw_interface.hh
@@ -1,217 +1,237 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef FFTW_INTERFACE
#define FFTW_INTERFACE
/* -------------------------------------------------------------------------- */
#include
+#include
+#include
#include
#include
#ifdef USE_MPI
#include
#endif
namespace fftw {
template
struct helper;
template <>
struct helper {
using complex = fftw_complex;
using plan = fftw_plan;
static auto alloc_real(std::size_t size) { return fftw_alloc_real(size); }
static auto alloc_complex(std::size_t size) {
return fftw_alloc_complex(size);
}
};
template <>
struct helper {
using complex = fftwl_complex;
using plan = fftwl_plan;
static auto alloc_real(std::size_t size) { return fftwl_alloc_real(size); }
static auto alloc_complex(std::size_t size) {
return fftwl_alloc_complex(size);
}
};
template
inline auto free(T* ptr) {
fftw_free(ptr);
}
inline auto destroy(fftw_plan plan) { fftw_destroy_plan(plan); }
inline auto destroy(fftwl_plan plan) { fftwl_destroy_plan(plan); }
inline auto init_threads() { return fftw_init_threads(); }
inline auto plan_with_nthreads(int nthreads) {
return fftw_plan_with_nthreads(nthreads);
}
inline auto cleanup_threads() { return fftw_cleanup_threads(); }
-inline auto mpi_init() {
+
+namespace mpi {
+inline auto init() {
#ifdef USE_MPI
return fftw_mpi_init();
#endif
}
-inline auto mpi_cleanup() {
+inline auto cleanup() {
#ifdef USE_MPI
return fftw_mpi_cleanup();
#endif
}
+inline auto local_size_many(int rank, const std::ptrdiff_t* size,
+ std::ptrdiff_t howmany) {
+#ifdef USE_MPI
+ std::ptrdiff_t local_n0, local_n0_offset;
+ auto res =
+ fftw_mpi_local_size_many(rank, size, howmany, FFTW_MPI_DEFAULT_BLOCK,
+ MPI_COMM_WORLD, &local_n0, &local_n0_offset);
+ return std::make_tuple(res, local_n0, local_n0_offset);
+#else
+ return std::make_tuple(howmany * std::accumulate(size, size + rank,
+ std::ptrdiff_t{1},
+ std::multiplies()),
+ size[0], 0);
+#endif
+}
+} // namespace mpi
/// Holder type for fftw plans
template
struct plan {
typename helper::plan _plan;
/// Create from plan
explicit plan(typename helper::plan _plan = nullptr) : _plan(_plan) {}
/// Move constructor to avoid accidental plan destruction
plan(plan&& o) : _plan(std::exchange(o._plan, nullptr)) {}
/// Move operator
plan& operator=(plan&& o) {
_plan = std::exchange(o._plan, nullptr);
return *this;
}
/// Destroy plan
~plan() noexcept {
if (_plan)
fftw::destroy(_plan);
}
/// For seamless use with fftw api
operator typename helper::plan() const { return _plan; }
};
/// RAII helper for fftw_free
template
struct ptr {
T* _ptr;
~ptr() noexcept {
if (_ptr)
fftw::free(_ptr);
}
operator T*() { return _ptr; }
};
/* -------------------------------------------------------------------------- */
inline auto plan_many_forward(int rank, const int* n, int howmany, double* in,
const int* inembed, int istride, int idist,
fftw_complex* out, const int* onembed,
int ostride, int odist, unsigned flags) {
return fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist,
out, onembed, ostride, odist, flags);
}
inline auto plan_many_backward(int rank, const int* n, int howmany,
fftw_complex* in, const int* inembed,
int istride, int idist, double* out,
const int* onembed, int ostride, int odist,
unsigned flags) {
return fftw_plan_many_dft_c2r(rank, n, howmany, in, inembed, istride, idist,
out, onembed, ostride, odist, flags);
}
inline auto plan_1d_forward(int n, double* in, fftw_complex* out,
unsigned flags) {
return fftw_plan_dft_r2c_1d(n, in, out, flags);
}
inline auto plan_1d_backward(int n, fftw_complex* in, double* out,
unsigned flags) {
return fftw_plan_dft_c2r_1d(n, in, out, flags);
}
inline auto plan_2d_forward(int n0, int n1, double* in, fftw_complex* out,
unsigned flags) {
return fftw_plan_dft_r2c_2d(n0, n1, in, out, flags);
}
inline auto plan_2d_backward(int n0, int n1, fftw_complex* out, double* in,
unsigned flags) {
return fftw_plan_dft_c2r_2d(n0, n1, out, in, flags);
}
inline auto execute(fftw_plan plan) { fftw_execute(plan); }
inline auto execute(fftw_plan plan, double* in, fftw_complex* out) {
fftw_execute_dft_r2c(plan, in, out);
}
inline auto execute(fftw_plan plan, fftw_complex* in, double* out) {
fftw_execute_dft_c2r(plan, in, out);
}
/* -------------------------------------------------------------------------- */
inline auto plan_many_forward(int rank, const int* n, int howmany,
long double* in, const int* inembed, int istride,
int idist, fftwl_complex* out, const int* onembed,
int ostride, int odist, unsigned flags) {
return fftwl_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist,
out, onembed, ostride, odist, flags);
}
inline auto plan_many_backward(int rank, const int* n, int howmany,
fftwl_complex* in, const int* inembed,
int istride, int idist, long double* out,
const int* onembed, int ostride, int odist,
unsigned flags) {
return fftwl_plan_many_dft_c2r(rank, n, howmany, in, inembed, istride, idist,
out, onembed, ostride, odist, flags);
}
inline auto plan_1d_forward(int n, long double* in, fftwl_complex* out,
unsigned flags) {
return fftwl_plan_dft_r2c_1d(n, in, out, flags);
}
inline auto plan_1d_backward(int n, fftwl_complex* in, long double* out,
unsigned flags) {
return fftwl_plan_dft_c2r_1d(n, in, out, flags);
}
inline auto plan_2d_forward(int n0, int n1, long double* in, fftwl_complex* out,
unsigned flags) {
return fftwl_plan_dft_r2c_2d(n0, n1, in, out, flags);
}
inline auto plan_2d_backward(int n0, int n1, fftwl_complex* out,
long double* in, unsigned flags) {
return fftwl_plan_dft_c2r_2d(n0, n1, out, in, flags);
}
inline auto execute(fftwl_plan plan) { fftwl_execute(plan); }
inline auto execute(fftwl_plan plan, long double* in, fftwl_complex* out) {
fftwl_execute_dft_r2c(plan, in, out);
}
inline auto execute(fftwl_plan plan, fftwl_complex* in, long double* out) {
fftwl_execute_dft_c2r(plan, in, out);
}
} // namespace fftw
#endif // FFTW_INTERFACE
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index fd3b74f..f7bd770 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,338 +1,339 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __GRID_BASE_HH__
-#define __GRID_BASE_HH__
+#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
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Dimension agnostic grid with components stored per points
*/
template
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)) {}
/// Destructor
virtual ~GridBase() = default;
/* ------------------------------------------------------------------------ */
/* Iterator class */
/* ------------------------------------------------------------------------ */
public:
using iterator = iterator_::iterator;
using const_iterator = iterator_::iterator;
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(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& 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 \
GridBase& operator op(const GridBase& other)
VEC_OPERATOR(+=);
VEC_OPERATOR(*=);
VEC_OPERATOR(-=);
VEC_OPERATOR(/=);
#undef VEC_OPERATOR
template
T dot(const GridBase& 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 \
GridBase& operator op(const StaticArray& 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(globalDataSize());
- }
+ T mean() const { return sum() / static_cast(globalDataSize()); }
T var() const;
/* ------------------------------------------------------------------------ */
/* 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
void copy(const GridBase& other) {
if (other.dataSize() != this->dataSize())
this->resize(other.dataSize());
thrust::copy(other.begin(), other.end(), this->begin());
nb_components = other.nb_components;
}
template
void move(GridBase&& 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 data;
UInt nb_components = 1;
};
/* -------------------------------------------------------------------------- */
template
void GridBase::uniformSetComponents(const GridBase& 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 \
inline GridBase& GridBase::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 \
template \
inline GridBase& GridBase::operator op(const GridBase& 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 \
template \
GridBase& GridBase::operator op(const StaticArray& b) { \
TAMAAS_ASSERT(this->dataSize() % size == 0, \
"Broadcast operator cannot broadcast" \
<< size << " on array of size " << this->dataSize()); \
Tensor::type, size> a(b); \
Loop::loop([&a] CUDA_LAMBDA(VectorProxy g) { g op a; }, \
range>(*this)); \
return *this; \
}
BROADCAST_OPERATOR_IMPL(+=)
BROADCAST_OPERATOR_IMPL(-=)
BROADCAST_OPERATOR_IMPL(*=)
BROADCAST_OPERATOR_IMPL(/=)
BROADCAST_OPERATOR_IMPL(=)
#undef BROADCAST_OPERATOR_IMPL
/* -------------------------------------------------------------------------- */
template
inline T GridBase::min() const {
// const auto id = [] CUDA_LAMBDA (const T&x){return x;};
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::max() const {
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::sum() const {
return Loop::reduce([] CUDA_LAMBDA(const T& x) { return x; },
*this);
}
/* -------------------------------------------------------------------------- */
template
inline T GridBase::var() const {
const T mu = mean();
const T var = Loop::reduce(
[mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
- std::cout << "var = " << var << '\n';
return var / (globalDataSize() - 1);
}
/* -------------------------------------------------------------------------- */
template
template
T GridBase::dot(const GridBase& other) const {
return Loop::reduce(
[] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
}
} // namespace tamaas
#endif // __GRID_BASE_HH__
diff --git a/src/core/tamaas.cpp b/src/core/tamaas.cpp
index 3ec23cb..e1e802c 100644
--- a/src/core/tamaas.cpp
+++ b/src/core/tamaas.cpp
@@ -1,78 +1,78 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "tamaas.hh"
#include "fftw_interface.hh"
#include "mpi_interface.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
void initialize(UInt num_threads) {
mpi::thread provided = mpi::thread::single;
if (not mpi::initialized())
mpi::init_thread(nullptr, nullptr, mpi::thread::single, &provided);
bool should_init_threads = (provided > mpi::thread::single);
if (num_threads) {
omp_set_num_threads(num_threads); // set user-defined number of threads
}
if (should_init_threads and (not fftw::init_threads())) {
TAMAAS_EXCEPTION("FFTW could not initialize threads!");
}
if (mpi::initialized())
- fftw::mpi_init();
+ fftw::mpi::init();
if (should_init_threads) {
fftw::plan_with_nthreads(omp_get_max_threads());
}
}
/* -------------------------------------------------------------------------- */
void finalize() {
if (not mpi::finalized()) {
fftw::cleanup_threads();
- fftw::mpi_cleanup();
- MPI_Finalize();
+ fftw::mpi::cleanup();
+ mpi::finalize();
}
}
namespace {
/// Manager for initialize + finalize
struct entry_exit_points {
entry_exit_points() { initialize(); }
~entry_exit_points() { finalize(); }
static const entry_exit_points singleton;
};
const entry_exit_points entry_exit_points::singleton;
} // namespace
} // namespace tamaas
diff --git a/src/core/tamaas.cpp b/src/mpi/partitioner.hh
similarity index 55%
copy from src/core/tamaas.cpp
copy to src/mpi/partitioner.hh
index 3ec23cb..17a9c8d 100644
--- a/src/core/tamaas.cpp
+++ b/src/mpi/partitioner.hh
@@ -1,78 +1,63 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
-
-#include "tamaas.hh"
+#ifndef PARTITIONER_HH
+#define PARTITIONER_HH
+/* -------------------------------------------------------------------------- */
#include "fftw_interface.hh"
#include "mpi_interface.hh"
-#include
-
+#include "tamaas.hh"
+#include
+#include
+#include
/* -------------------------------------------------------------------------- */
-
namespace tamaas {
-void initialize(UInt num_threads) {
- mpi::thread provided = mpi::thread::single;
-
- if (not mpi::initialized())
- mpi::init_thread(nullptr, nullptr, mpi::thread::single, &provided);
-
- bool should_init_threads = (provided > mpi::thread::single);
-
- if (num_threads) {
- omp_set_num_threads(num_threads); // set user-defined number of threads
+template
+struct Partitioner {
+ static decltype(auto) global_size(std::array local) {
+ local.front() = mpi::allreduce(local.front());
+ return local;
}
- if (should_init_threads and (not fftw::init_threads())) {
- TAMAAS_EXCEPTION("FFTW could not initialize threads!");
+ static decltype(auto) local_size(std::array global) {
+ auto tup =
+ fftw::mpi::local_size_many(dim, cast_size(global).data(), 1);
+ global.front() = static_cast(std::get<1>(tup));
+ return global;
}
- if (mpi::initialized())
- fftw::mpi_init();
-
- if (should_init_threads) {
- fftw::plan_with_nthreads(omp_get_max_threads());
+ static decltype(auto) local_offset(const std::array& global) {
+ auto tup =
+ fftw::mpi::local_size_many(dim, cast_size(global).data(), 1);
+ return static_cast(std::get<2>(tup));
}
-}
-/* -------------------------------------------------------------------------- */
-
-void finalize() {
- if (not mpi::finalized()) {
- fftw::cleanup_threads();
- fftw::mpi_cleanup();
- MPI_Finalize();
+ static decltype(auto) cast_size(const std::array& s) {
+ std::array n;
+ std::copy(s.begin(), s.end(), n.begin());
+ return n;
}
-}
-
-namespace {
-/// Manager for initialize + finalize
-struct entry_exit_points {
- entry_exit_points() { initialize(); }
- ~entry_exit_points() { finalize(); }
- static const entry_exit_points singleton;
};
-const entry_exit_points entry_exit_points::singleton;
-} // namespace
-
} // namespace tamaas
+#endif
diff --git a/src/surface/surface_generator.cpp b/src/surface/surface_generator.cpp
index 5bc3485..ea96b63 100644
--- a/src/surface/surface_generator.cpp
+++ b/src/surface/surface_generator.cpp
@@ -1,49 +1,60 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "surface_generator.hh"
+#include "partitioner.hh"
+#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
-void SurfaceGenerator::setSizes(const std::array& n) {
- grid.resize(n);
+void SurfaceGenerator::setSizes(std::array n) {
+ global_size = n;
+ grid.resize(Partitioner::local_size(std::move(n)));
}
/* -------------------------------------------------------------------------- */
template
void SurfaceGenerator::setSizes(const UInt n[dim]) {
std::array size;
- std::copy(n, n + dim, size.begin());
- setSizes(size);
+ std::copy_n(n, dim, size.begin());
+ setSizes(std::move(size));
+}
+
+/* -------------------------------------------------------------------------- */
+
+template
+void SurfaceGenerator::setRandomSeed(long seed) {
+ random_seed = seed;
+ random_seed += Partitioner::local_offset(global_size);
}
/* -------------------------------------------------------------------------- */
template class SurfaceGenerator<1>;
template class SurfaceGenerator<2>;
} // namespace tamaas
diff --git a/src/surface/surface_generator.hh b/src/surface/surface_generator.hh
index 830e5df..3dc8a48 100644
--- a/src/surface/surface_generator.hh
+++ b/src/surface/surface_generator.hh
@@ -1,63 +1,66 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __SURFACE_GENERATOR_HH__
-#define __SURFACE_GENERATOR_HH__
+#ifndef SURFACE_GENERATOR_HH
+#define SURFACE_GENERATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
+#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Class generating random surfaces
template
class SurfaceGenerator {
public:
/// Default constructor
SurfaceGenerator() = default;
/// Default destructor
virtual ~SurfaceGenerator() = default;
public:
/// Build surface profile (array of heights)
virtual Grid& buildSurface() = 0;
/// Set surface sizes
- void setSizes(const std::array& n);
+ void setSizes(std::array n);
/// Set surface sizes
void setSizes(const UInt n[dim]);
- TAMAAS_ACCESSOR(random_seed, long, RandomSeed);
+ long getRandomSeed() const { return random_seed; }
+ void setRandomSeed(long seed);
/// Maintaining old interface
friend class SurfaceGeneratorFilterFFT;
protected:
Grid grid;
+ std::array global_size;
long random_seed = 0;
};
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // __SURFACE_GENERATOR_HH__
diff --git a/src/surface/surface_generator_filter.cpp b/src/surface/surface_generator_filter.cpp
index ab3496e..b5767e6 100644
--- a/src/surface/surface_generator_filter.cpp
+++ b/src/surface/surface_generator_filter.cpp
@@ -1,73 +1,73 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "surface_generator_filter.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
void SurfaceGeneratorFilter::applyFilterOnSource() {
GridHermitian fft_white_noise(filter_coefficients.sizes(), 1);
engine.forward(this->white_noise, fft_white_noise);
fft_white_noise *= filter_coefficients;
engine.backward(this->grid, fft_white_noise);
}
/* -------------------------------------------------------------------------- */
template
Grid& SurfaceGeneratorFilter::buildSurface() {
if (this->grid.dataSize() == 0)
TAMAAS_EXCEPTION("the size of the grid is zero, did you call setSizes() ?");
if (this->filter == nullptr)
TAMAAS_EXCEPTION("filter is null, did you call setFilter() ?");
// Resizing arrays
this->white_noise.resize(this->grid.sizes());
auto hermitian_dim =
GridHermitian::hermitianDimensions(this->grid.sizes());
this->filter_coefficients.resize(hermitian_dim);
/// Generating white noise with gaussian distribution
this->generateWhiteNoise>();
/// Applying filters
filter->computeFilter(this->filter_coefficients);
this->applyFilterOnSource();
/// Normalizing surface for hrms independent of number of points
- this->grid *= std::sqrt(this->grid.dataSize());
+ this->grid *= std::sqrt(this->grid.globalDataSize());
return this->grid;
}
/* -------------------------------------------------------------------------- */
template class SurfaceGeneratorFilter<1>;
template class SurfaceGeneratorFilter<2>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/surface/surface_generator_filter_fft.hh b/src/surface/surface_generator_filter_fft.hh
index d730ddd..0c50f60 100644
--- a/src/surface/surface_generator_filter_fft.hh
+++ b/src/surface/surface_generator_filter_fft.hh
@@ -1,74 +1,74 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef SURFACE_GENERATOR_FILTER_FFT_H
#define SURFACE_GENERATOR_FILTER_FFT_H
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include "surface_generator_filter.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Legacy class to maintain interface
class SurfaceGeneratorFilterFFT {
public:
//! can construct a generator with default values
SurfaceGeneratorFilterFFT();
public:
//! get high frequency cut
UInt& getQ2() { return pl_filter->getQ2(); }
//! get low frequency cut
UInt& getQ0() { return pl_filter->getQ0(); }
//! get roll frequency cut
UInt& getQ1() { return pl_filter->getQ1(); }
//! get hurst
Real& getHurst() { return pl_filter->getHurst(); }
/// compute analytical RMS
Real analyticalRMS() const { return pl_filter->rmsHeights(); }
/// Swig does not wrap the numpy otherwise
Grid& buildSurface();
/// Get gridsize
UInt& getGridSize() { return grid_size; }
/// Get random seed
- long& getRandomSeed() { return generator.getRandomSeed(); }
+ long& getRandomSeed() { return generator.random_seed; }
/// Get RMS
Real& getRMS() { return rms; }
/// Contrain RMS
Real constrainRMS();
/// Init
void Init() {}
protected:
SurfaceGeneratorFilter<2> generator;
std::shared_ptr> pl_filter = nullptr;
UInt grid_size;
Real rms;
};
} // namespace tamaas
#endif
diff --git a/src/surface/surface_generator_random_phase.cpp b/src/surface/surface_generator_random_phase.cpp
index b9af274..fbdc690 100644
--- a/src/surface/surface_generator_random_phase.cpp
+++ b/src/surface/surface_generator_random_phase.cpp
@@ -1,67 +1,67 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "surface_generator_random_phase.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
Grid& SurfaceGeneratorRandomPhase::buildSurface() {
if (this->grid.dataSize() == 0)
TAMAAS_EXCEPTION("the size of the grid is zero, did you call setSizes() ?");
if (this->filter == nullptr)
- TAMAAS_EXCEPTION("filter is null, did you call setFilter() ?");
+ TAMAAS_EXCEPTION("spectrum is null, did you call setSpectrum() ?");
// Resizing arrays
auto hermitian_dim =
GridHermitian::hermitianDimensions(this->grid.sizes());
this->white_noise.resize(hermitian_dim);
this->filter_coefficients.resize(hermitian_dim);
/// Generating white noise with uniform distribution for phase
this->template generateWhiteNoise>();
/// Applying filters
this->filter->computeFilter(this->filter_coefficients);
Loop::loop([](Complex& coeff,
Real& phase) { coeff *= thrust::polar(1., 2 * M_PI * phase); },
this->filter_coefficients, this->white_noise);
this->engine.backward(this->grid, this->filter_coefficients);
/// Normalizing surface for hrms independent of number of points
- this->grid *= this->grid.dataSize();
+ this->grid *= this->grid.globalDataSize();
return this->grid;
}
/* -------------------------------------------------------------------------- */
template class SurfaceGeneratorRandomPhase<1>;
template class SurfaceGeneratorRandomPhase<2>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas