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