diff --git a/.vscode/extensions.json b/.vscode/extensions.json index ddcd3e2..4af38e8 100644 --- a/.vscode/extensions.json +++ b/.vscode/extensions.json @@ -1,7 +1,6 @@ { "recommendations": [ "ms-vscode.cpptools", - "twxs.cmake", "ms-vscode.cmake-tools", "ms-python.python", "codecreator.indent-to-bracket", diff --git a/.vscode/settings.json b/.vscode/settings.json index 6b04c16..b8347a0 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -106,5 +106,5 @@ "__bits": "cpp", "__verbose_abort": "cpp" }, - "cmake.configureOnOpen": true + "cmake.configureOnOpen": false } diff --git a/include/flyft/cartesian_mesh.h b/include/flyft/cartesian_mesh.h index ee8fa06..b497e35 100644 --- a/include/flyft/cartesian_mesh.h +++ b/include/flyft/cartesian_mesh.h @@ -6,20 +6,43 @@ namespace flyft { +//! Mesh with Cartesian position coordinates. class CartesianMesh : public Mesh { public: - CartesianMesh(double lower_bound, + //! Constructor. + /*! + * \param num_orientation_dim Number of orientation coordinates. + * \param shape Number of cells along each dimension. + * \param lower_bound Lower bound of the mesh. + * \param upper_bound Upper bound of the mesh. + * \param lower_bc Boundary condition at position coordinate lower bound. + * \param upper_bc Boundary condition at position coordinate upper bound. + * \param area Cross-sectional area, tangent to position coordinate. + */ + CartesianMesh(int num_orientation_dim, + int* shape, + double lower_bound, double upper_bound, - int shape, BoundaryType lower_bc, BoundaryType upper_bc, double area); - double area(int i) const override; - double volume() const override; - double volume(int i) const override; - double gradient(int idx, double f_lo, double f_hi) const override; + //! Copy constructor. + CartesianMesh(const CartesianMesh& other); + + //! Copy assignment operator. + CartesianMesh& operator=(const CartesianMesh& other); + + //! Move constructor. + CartesianMesh(const CartesianMesh&& other); + + //! Move assignment operator. + CartesianMesh& operator=(const CartesianMesh&& other); + + // double area(int i) const override; + double cell_position_volume(const int* cell) const override; + // double gradient(int idx, double f_lo, double f_hi) const override; protected: std::shared_ptr clone() const override; diff --git a/include/flyft/complex_field.h b/include/flyft/complex_field.h new file mode 100644 index 0000000..e217551 --- /dev/null +++ b/include/flyft/complex_field.h @@ -0,0 +1,20 @@ +#ifndef FLYFT_COMPLEX_FIELD_H_ +#define FLYFT_COMPLEX_FIELD_H_ + +#include "flyft/data_array.h" + +#include + +namespace flyft + { + +//! Complex-valued field. +/*! + * A ComplexField represents a complex-valued quantity discretized onto a coordinate space. It is a + * multidimensional DataArray and requires another object to define the coordinates. + */ +using ComplexField = DataArray, 4>; + + } // namespace flyft + +#endif // FLYFT_COMPLEX_FIELD_H_ diff --git a/include/flyft/data_array.h b/include/flyft/data_array.h new file mode 100644 index 0000000..d671e86 --- /dev/null +++ b/include/flyft/data_array.h @@ -0,0 +1,346 @@ +#ifndef FLYFT_DATA_ARRAY_H_ +#define FLYFT_DATA_ARRAY_H_ + +#include "flyft/data_layout.h" +#include "flyft/data_view.h" +#include "flyft/tracked_object.h" + +#include + +namespace flyft + { + +//! Multidimensional array. +/*! + * A DataArray is a multidimensional array. It allocates and stores its own data consistent with a + * certain number of dimensions and shape, as defined by a DataLayout. + * + * The \b full shape of the DataArray includes two sub-shapes: the shape and the buffer shape. The + * full shape is the shape plus one buffer shape on both sides. The intention is that the shape be + * used for locally owned data and the buffers be used for data owned by other processors when + * arrays are broken across processors in parallel simulations. + * + * A new buffer shape can be requested, and the buffer shape along each dimension is the largest + * value requested. + * + * A DataArray can be reshaped. If the buffer shape is changed but the shape is not, the data stored + * in the shape will be copied. + * + * A DataArray is a TrackedObject. This means that potential modifications to its values are tracked + * when a view of the data is taken. For performance elsewhere, it may be important to use + * const_view when the data is only being read. + * + * The maximum dimension that can be represented by the layout is set by the template parameter + * \a N. + */ +template +class DataArray : public TrackedObject + { + public: + using View = DataView; + using ConstantView = DataView; + using Iterator = typename View::Iterator; + using ConstantIterator = typename ConstantView::Iterator; + + // no default constructor + DataArray() = delete; + + //! Constructor. + /*! + * \param num_dimensions Number of dimensions. + * \param shape Number of elements per dimension. + */ + DataArray(int num_dimensions, const int* shape); + + //! Constructor with buffer. + /*! + * \param num_dimensions Number of dimensions. + * \param shape Number of elements per dimension. + * \param buffer_shape Number of elements in buffers per dimension. + */ + DataArray(int num_dimensions, const int* shape, const int* buffer_shape); + + // noncopyable / nonmovable + DataArray(const DataArray&) = delete; + DataArray(DataArray&&) = delete; + DataArray& operator=(const DataArray&) = delete; + DataArray& operator=(DataArray&&) = delete; + + //! Destructor. + ~DataArray(); + + //! Access an element. + /*! + * \param multi_index Multidimensional index. + * \return Element. + * + * Accessing this way allows the element to be modified, so it will advance the tracked state + * of the array. + */ + T& operator()(const int* multi_index); + + //! Access an element. + /*! + * \param multi_index Multidimensional index. + * \return Element. + */ + const T& operator()(const int* multi_index) const; + + //! Number of dimensions. + int num_dimensions() const; + + //! Number of elements per dimension. + const int* shape() const; + + //! Number of elements in each buffer per dimension. + const int* buffer_shape() const; + + //! Total number of elements per dimension. + /*! + * If there are $S$ elements in the shape and $B$ elements in the buffer for a dimension, then + * the total number of elements for the dimension is $F = S + 2B$. + */ + const int* full_shape() const; + + //! View of elements in the shape. + /*! + * This view allows the data array to be modified, so the tracked state of the array will be + * advanced when it is created. + */ + View view(); + + //! Read-only view of elements in the shape. + ConstantView const_view() const; + + //! View of elements in the full shape. + /*! + * This view allows the data array to be modified, so the tracked state of the array will be + * advanced when it is created. + */ + View full_view(); + + //! Read-only view of elements in the full shape. + ConstantView const_full_view() const; + + //! Reshape the data array. + /*! + * \param num_dimensions Number of dimensions. + * \param shape Number of elements per dimension. + * \param buffer_shape Number of buffer elements per dimension. + * + * If the shape stays the same but the buffer shape changes, the data in the shape will be + * copied to the new array. Otherwise, all data will be destroyed. + */ + void reshape(int num_dimensions, const int* shape, const int* buffer_shape); + + //! Force a new buffer shape. + /*! + * \param buffer_shape Number of buffer elements per dimension. + * + * The buffer shape is immediately resized. + */ + void setBuffer(const int* buffer_shape); + + //! Request a new buffer shape. + /*! + * \param buffer_shape Number of buffer elements per dimension. + * + * The buffer shape only changes if it increases in at least one dimension. + */ + void requestBuffer(const int* buffer_shape); + + private: + T* data_; //!< Data. + int shape_[N]; //!< Number of elements per dimension. + int buffer_shape_[N]; //!< Number of elements in buffer per dimension. + typename View::Layout layout_; //!< Layout of the data array. + }; + +template +DataArray::DataArray(int num_dimensions, const int* shape) + : DataArray(num_dimensions, shape, nullptr) + { + } + +template +DataArray::DataArray(int num_dimensions, const int* shape, const int* buffer_shape) + : shape_ {0}, buffer_shape_ {0} + { + assert(num_dimensions > 0 && num_dimensions <= N); + assert(shape); + reshape(num_dimensions, shape, buffer_shape); + } + +template +DataArray::~DataArray() + { + free(data_); + data_ = nullptr; + } + +template +T& DataArray::operator()(const int* multi_index) + { + assert(multi_index); + return view()(multi_index); + } + +template +const T& DataArray::operator()(const int* multi_index) const + { + assert(multi_index); + return const_view()(multi_index); + } + +template +int DataArray::num_dimensions() const + { + return layout_.num_dimensions(); + } + +template +const int* DataArray::shape() const + { + return shape_; + } + +template +const int* DataArray::buffer_shape() const + { + return buffer_shape_; + } + +template +const int* DataArray::full_shape() const + { + return layout_.shape(); + } + +template +typename DataArray::View DataArray::view() + { + token_.stageAndCommit(); + return View(data_, layout_, buffer_shape_, shape_); + } + +template +typename DataArray::ConstantView DataArray::const_view() const + { + return ConstantView(static_cast(data_), layout_, buffer_shape_, shape_); + } + +template +typename DataArray::View DataArray::full_view() + { + token_.stageAndCommit(); + return View(data_, layout_); + } + +template +typename DataArray::ConstantView DataArray::const_full_view() const + { + return ConstantView(data_, layout_); + } + +template +void DataArray::reshape(int num_dimensions, const int* shape, const int* buffer_shape) + { + assert(num_dimensions > 0 && num_dimensions <= N); + + bool same_shape = true; + bool same_buffer_shape = true; + std::vector full_shape(num_dimensions); + for (int dim = 0; dim < num_dimensions; ++dim) + { + if (shape && shape[dim] != shape_[dim]) + { + same_shape = false; + shape_[dim] = shape[dim]; + assert(shape_[dim] > 0); + } + + if (buffer_shape && buffer_shape[dim] != buffer_shape_[dim]) + { + same_buffer_shape = false; + buffer_shape_[dim] = buffer_shape[dim]; + assert(buffer_shape_[dim] >= 0); + } + + full_shape[dim] = shape_[dim] + 2 * buffer_shape_[dim]; + } + + if (!(same_shape && same_buffer_shape)) + { + typename View::Layout layout(num_dimensions, full_shape.data()); + + // if data is alloc'd, decide what to do with it + if (data_) + { + if (same_shape) + { + // data shape is the same but buffer changed, copy + T* tmp = static_cast(malloc(layout.size() * sizeof(T))); + std::fill(tmp, tmp + layout.size(), T(0)); + + // copy the data from old to new + auto view_old = const_view(); + View view_new(tmp, layout, buffer_shape, shape); + std::copy(view_old.begin(), view_old.end(), view_new.begin()); + std::swap(data_, tmp); + free(tmp); + } + else if (layout.size() == layout_.size() && layout.size() > 0) + { + // total shape is still the same, just reset the data + std::fill(data_, data_ + layout.size(), T(0)); + } + else + { + // shape is different, wipe it out and realloc + free(data_); + data_ = nullptr; + } + } + + // if data is still not alloc'd, make it so + if (!data_ && layout.size() > 0) + { + data_ = static_cast(malloc(layout.size() * sizeof(T))); + std::fill(data_, data_ + layout.size(), T(0)); + } + layout_ = layout; + token_.stageAndCommit(); + } + } + +template +void DataArray::setBuffer(const int* buffer_shape) + { + reshape(layout_.num_dimensions(), shape_, buffer_shape); + } + +template +void DataArray::requestBuffer(const int* buffer_shape) + { + // don't do anything if buffer shape is null + if (!buffer_shape) + return; + + bool any_increase = false; + for (int dim = 0; !any_increase && dim < layout_.num_dimensions(); ++dim) + { + if (buffer_shape[dim] > buffer_shape_[dim]) + { + any_increase = true; + } + } + + if (any_increase) + { + setBuffer(buffer_shape); + } + } + + } // namespace flyft + +#endif // FLYFT_DATA_ARRAY_H_ diff --git a/include/flyft/data_layout.h b/include/flyft/data_layout.h index ca04fb9..c424c67 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -1,28 +1,250 @@ #ifndef FLYFT_DATA_LAYOUT_H_ #define FLYFT_DATA_LAYOUT_H_ -#include -#include +#include +#include namespace flyft { + +//! Memory layout of a multidimensional array. +/*! + * A multidimensional array has a number of dimensions \a n. Its shape is the number of elements + * along each dimension, and its size is the total number of elements (product of the shape). An + * element of the array can be accessed by a multidimensional index, which is a length \a n array. + * The array is layed out contiguously in one-dimensional memory using generalized row-major + * ordering, meaning that the first index is the "slowest" varying and the last index is the + * "fastest" varying. + * + * The maximum dimension that can be represented by the layout is set by the template parameter + * \a N. + */ +template class DataLayout { public: + //! Empty constructor DataLayout(); - explicit DataLayout(int shape_); - int operator()(int idx) const; - int shape() const; - int size() const; + //! Constructor + /*! + * \param num_dimensions Number of dimensions. + * \param shape Number of elements per dimension. + */ + DataLayout(int num_dimensions, const int* shape); + + //! Map a multidimensional index to a one-dimensional index. + /*! + * \param multi_index Multidimensional index. + * \return One-dimensional index. + */ + size_t operator()(const int* multi_index) const; + + //! Map a multidimensional index offset by another into a one-dimensional index. + /*! + * \param multi_index Multidimensional index. + * \param offset Offset for multidimensional index. + * \return One-dimensional index. + * + * The total multidimensional index is assumed to be valid. + */ + size_t operator()(const int* multi_index, const int* offset) const; + + //! Map a one-dimensional index to a multidimensional index. + /*! + * \param multi_index Multidimensional index. + * \param flat_index One-dimensional index. + */ + void operator()(int* multi_index, size_t flat_index) const; + + //! Map a one-dimensional index to a multidimensional index relative to an offset. + /*! + * \param multi_index Multidimensional index. + * \param flat_index One-dimensional index. + * \param offset Offset for multidimensional index. + */ + void operator()(int* multi_index, size_t flat_index, const int* offset) const; + + //! Number of dimensions. + int num_dimensions() const; + + //! Number of elements per dimension. + const int* shape() const; + + //! Total number of elements. + size_t size() const; + + //! Test if two layouts are equal. + template + bool operator==(const DataLayout& other) const; + + //! Test if two layouts are not equal. + template + bool operator!=(const DataLayout& other) const; - bool operator==(const DataLayout& other) const; - bool operator!=(const DataLayout& other) const; + static int constexpr MAX_NUM_DIMENSIONS = N; //!< Maximum number of dimensions. private: - int shape_; + int num_dimensions_; //!< Number of dimensions. + int shape_[N]; //!< Number of elements per dimension. + size_t stride_[N]; //!< Number of elements to advance when incrementing index. + size_t size_; //!< Total number of elements. }; +template +DataLayout::DataLayout() : num_dimensions_(0), shape_ {0}, stride_ {0}, size_(0) + { + } + +template +DataLayout::DataLayout(int num_dimensions, const int* shape) : num_dimensions_(num_dimensions) + { + assert(num_dimensions > 0 && num_dimensions_ <= N); + assert(shape); + + // row-major ordering + size_t stride = 1; + for (int dim = num_dimensions - 1; dim >= 0; --dim) + { + assert(shape[dim] >= 1); + + shape_[dim] = shape[dim]; + stride_[dim] = stride; + stride *= shape_[dim]; + } + + // final value of stride will be the full size + size_ = stride; + } + +template +size_t DataLayout::operator()(const int* multi_index) const + { + assert(multi_index); + + size_t flat_index; + if (num_dimensions_ > 1) + { + flat_index = 0; + for (int dim = 0; dim < num_dimensions_; ++dim) + { + flat_index += multi_index[dim] * stride_[dim]; + } + } + else + { + flat_index = multi_index[0]; + } + return flat_index; + } + +template +size_t DataLayout::operator()(const int* multi_index, const int* offset) const + { + // use method without addition if there is no starting index. + if (!offset) + { + return this->operator()(multi_index); + } + + // otherwise, do calculation with starting offset added + assert(multi_index); + assert(offset); + size_t flat_index; + if (num_dimensions_ > 1) + { + flat_index = 0; + for (int dim = 0; dim < num_dimensions_; ++dim) + { + flat_index += (offset[dim] + multi_index[dim]) * stride_[dim]; + } + } + else + { + flat_index = offset[0] + multi_index[0]; + } + return flat_index; + } + +template +void DataLayout::operator()(int* multi_index, size_t flat_index) const + { + assert(multi_index); + + if (num_dimensions_ > 1) + { + // repeatedly floor divide the flat index by stride associated with it + size_t idx = flat_index; + for (int dim = 0; dim < num_dimensions_; ++dim) + { + multi_index[dim] = idx / stride_[dim]; + idx -= multi_index[dim] * stride_[dim]; + } + } + else + { + multi_index[0] = flat_index; + } + } + +template +void DataLayout::operator()(int* multi_index, size_t flat_index, const int* offset) const + { + this->operator()(multi_index, flat_index); + if (offset) + { + for (int dim = 0; dim < num_dimensions_; ++dim) + { + multi_index[dim] -= offset[dim]; + } + } + } + +template +int DataLayout::num_dimensions() const + { + return num_dimensions_; + } + +template +const int* DataLayout::shape() const + { + return shape_; + } + +template +size_t DataLayout::size() const + { + return size_; + } + +template +template +bool DataLayout::operator==(const DataLayout& other) const + { + // dimensionality & size must match as first pass + bool same_layout = (num_dimensions_ == other.num_dimensions_ && size_ == other.size_); + if (same_layout) + { + for (int dim = 0; dim < num_dimensions_ && same_layout; ++dim) + { + if (shape_[dim] != other.shape_[dim]) + { + same_layout = false; + } + } + } + + return same_layout; + } + +template +template +bool DataLayout::operator!=(const DataLayout& other) const + { + return !(*this == other); + } + } // namespace flyft #endif // FLYFT_DATA_LAYOUT_H_ diff --git a/include/flyft/data_view.h b/include/flyft/data_view.h index 6da3f61..a08eb54 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -3,144 +3,386 @@ #include "flyft/data_layout.h" +#include +#include +#include +#include + namespace flyft { -template +//! Accessor to a multidimensional array. +/*! + * A DataView interprets a pointer as a multidimensional array according to a particular DataLayout. + * It provides random access to elements by multidimensional index, as well as a forward iterator. + * + * The view can be scoped to a subset of only some of the valid multidimensional indexes. This is + * useful for accessing only the elements of the array that are of interest when the array is + * padded. When the view is scoped, negative indexes can be used and are interpreted relative to the + * starting multidimensional index for the scope. One use of negative indexes, in conjunction with + * padding, is to apply finite-difference stencils near edges without special-case handling. It is + * assumed that the resulting index is valid for the layout. + * + * The maximum dimension that can be represented by the layout is set by the template parameter + * \a N. + */ +template class DataView { public: - using value_type = typename std::remove_reference::type; - using pointer = value_type*; - using reference = value_type&; + using value_type = typename std::remove_reference::type; //; //!< Layout type class Iterator { public: - using iterator_category = std::bidirectional_iterator_tag; - using difference_type = int; - using value_type = DataView::value_type; - using pointer = DataView::pointer; - using reference = DataView::reference; + using iterator_category = std::forward_iterator_tag; //!< Iterator type + using value_type = DataView::value_type; //!< Data type + using difference_type = std::ptrdiff_t; //!< Difference type + using pointer = DataView::pointer; //!< Pointer to data type + using reference = DataView::reference; //!< Reference to data type - Iterator() : Iterator(DataView()) {} + //! Empty constructor + Iterator(); - explicit Iterator(const DataView& view) : Iterator(view, 0) {} + //! Construct from view. + /*! + * \param view Data view to iterate. + * + * This iterator starts at the zero multidimensional index. + */ + explicit Iterator(const DataView& view); - Iterator(const DataView& view, int current) : view_(view), current_(current) {} + //! Construct from view starting at multidimensional index. + /*! + * \param view Data view to iterate. + * \param current Multi-index to start iterator at. + */ + Iterator(const DataView& view, const int* current); - reference operator*() const - { - return view_(current_); - } + //! Dereference iterator. + reference operator*() const; - pointer operator->() const - { - return get(); - } + //! Access current data as pointer. + pointer operator->() const; - pointer get() const - { - return &view_(current_); - } + //! Get pointer to current data. + pointer get() const; - Iterator& operator++() - { - ++current_; - return *this; - } + //! Pre-increment. + /*! + * The last dimension of the multidimensional index is incremented by 1. + */ + Iterator& operator++(); - Iterator operator++(int) - { - Iterator tmp(*this); - ++current_; - return tmp; - } + //! Post-increment. + /*! + * The last dimension of the multidimensional index is incremented by 1. + */ + Iterator operator++(int); - Iterator& operator--() - { - --current_; - return *this; - } + //! In-place increment. + /*! + * \param n Number of elements to increment. + * + * The last dimension of the multidimensional index is incremented by n. + */ + Iterator& operator+=(difference_type n); - Iterator operator--(int) - { - Iterator tmp(*this); - --current_; - return tmp; - } + //! Move iterator forward. + /*! + * \param n Number of elements to increment. + * + * The last dimension of the multidimensional index is incremented by n. + */ + Iterator operator+(difference_type n); - bool operator==(const Iterator& other) const - { - return (get() == other.get()); - } + //! Check if two iterators are equivalent. + /*! + * \return True if both iterators point to the same data. + * + * To make the comparison faster, it is assumed both iterators were constructed with the + * same view. + */ + bool operator==(const Iterator& other) const; - bool operator!=(const Iterator& other) const - { - return !(*this == other); - } + //! Check if two iterators are not equivalent. + /*! + * \return True if both iterators don't point to the same data. + */ + bool operator!=(const Iterator& other) const; private: - DataView view_; - int current_; + DataView view_; //!< View being iterated over + int current_[N]; //!< Current multidimensional index }; - DataView() : DataView(nullptr, DataLayout()) {} + //! Empty constructor. + DataView(); - DataView(pointer data, const DataLayout& layout) : DataView(data, layout, 0, layout.shape()) {} + //! Constructor. + /*! + * \param data Pointer to the data. + * \param layout Layout of the data. + */ + DataView(pointer data, const Layout& layout); - DataView(pointer data, const DataLayout& layout, int start, int end) - : data_(data), layout_(layout), start_(start), end_(end) - { - } + //! Scoped-view constructor. + /*! + * \param data Pointer to the data. + * \param layout Layout of the data. + * \param start Lower bound of multidimensional indexes within scope. + * \param shape Shape of multidimensional indexes within scope, beginning with start. + */ + DataView(pointer data, const Layout& layout, const int* start, const int* shape); - reference operator()(int idx) const - { - return data_[layout_(start_ + idx)]; - } + //! Access an element of the data at a given multidimensional index. + /*! + * \param multi_index Multidimensional index to access. + * \return Data element. + */ + reference operator()(const int* multi_index) const; + + //! Number of dimensions. + int num_dimensions() const; + + //! Number of elements per dimension in the scope of the view. + const int* shape() const; + + //! Number of elements in the scope of the view. + size_t size() const; + + //! Test if view is not null + /*! + * \returns True if data pointer being viewed is not null. + */ + explicit operator bool() const; + + //! Make an iterator to the start of the data being viewed. + /*! + * \return Iterator to start of data. + */ + Iterator begin() const; + + //! Make and iterator to the end of the data being viewed. + /*! + * \return Iterator to end of data. + */ + Iterator end() const; + + //! Make an iterator to a particular multidimensional index. + /*! + * \return Iterator to specified multidimensional index. + */ + Iterator make_iterator(const int* multi_index) const; + + static int constexpr MAX_NUM_DIMENSIONS = N; //!< Maximum number of dimensions. - int shape() const + private: + pointer data_; //!< Data. + Layout layout_; //!< Multidimensional array layout. + int start_[N]; //!< Multi-index for start of scoped view. + int shape_[N]; //!< Number of elements per dimension of scoped view. + }; + +template +DataView::DataView() : data_(nullptr), layout_(Layout()), start_ {0}, shape_ {0} + { + } + +template +DataView::DataView(pointer data, const Layout& layout) : data_(data), layout_(layout) + { + assert(data); + assert(layout.num_dimensions() > 0 && layout.num_dimensions() <= N); + + const auto layout_shape = layout_.shape(); + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { - return end_ - start_; + start_[dim] = 0; + shape_[dim] = layout_shape[dim]; } + } - int size() const +template +DataView::DataView(pointer data, const Layout& layout, const int* start, const int* shape) + : data_(data), layout_(layout) + { + assert(data); + assert(layout.num_dimensions() > 0 && layout.num_dimensions() <= N); + assert(start); + assert(shape); + + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { - return shape(); + assert(start_[dim] >= 0); + assert(shape_[dim] > 0); + assert(start_[dim] + shape_[dim] <= layout.shape()[dim]); + start_[dim] = start[dim]; + shape_[dim] = shape[dim]; } + } - explicit operator bool() const +template +typename DataView::reference DataView::operator()(const int* multi_index) const + { + assert(multi_index); + return data_[layout_(multi_index, start_)]; + } + +template +int DataView::num_dimensions() const + { + return layout_.num_dimensions(); + } + +template +const int* DataView::shape() const + { + return shape_; + } + +template +size_t DataView::size() const + { + size_t size; + if (shape_) { - return (data_ != nullptr); + size = 1; + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) + { + size *= shape_[dim]; + } } - - Iterator begin() const + else { - return Iterator(*this, 0); + size = 0; } + return size; + } + +template +DataView::operator bool() const + { + return (data_ != nullptr); + } + +template +typename DataView::Iterator DataView::begin() const + { + return Iterator(*this); + } + +template +typename DataView::Iterator DataView::end() const + { + return Iterator(*this, shape_); + } + +template +typename DataView::Iterator DataView::make_iterator(const int* multi_index) const + { + assert(multi_index); + return Iterator(*this, multi_index); + } + +template +DataView::Iterator::Iterator() : view_(DataView()), current_ {0} + { + } + +template +DataView::Iterator::Iterator(const DataView& view) : view_(view), current_ {0} + { + assert(view); + } - Iterator end() const +template +DataView::Iterator::Iterator(const DataView& view, const int* current) : view_(view) + { + assert(view); + for (int dim = 0; dim < view_.num_dimensions(); ++dim) { - return Iterator(*this, shape()); + current_[dim] = (current) ? current[dim] : 0; } + } + +template +typename DataView::Iterator::reference DataView::Iterator::operator*() const + { + return view_(current_); + } + +template +typename DataView::Iterator::pointer DataView::Iterator::operator->() const + { + return get(); + } - bool operator==(const DataView& other) const +template +typename DataView::Iterator::pointer DataView::Iterator::get() const + { + return &view_(current_); + } + +template +typename DataView::Iterator& DataView::Iterator::operator++() + { + // increment the multi-index of the last dimension by 1 + int dim = view_.num_dimensions() - 1; + ++current_[dim]; + + // then, carry forward + const auto view_shape = view_.shape(); + while (current_[dim] == view_shape[dim] && dim >= 1) { - return (data_ == other.data_ && layout_ == other.layout_ && start_ == other.start_ - && end_ == other.end_); + current_[dim] = 0; + ++current_[dim--]; } - bool operator!=(const DataView& other) const + return *this; + } + +template +typename DataView::Iterator DataView::Iterator::operator++(int) + { + Iterator tmp(*this); + ++(*this); + return tmp; + } + +template +typename DataView::Iterator& DataView::Iterator::operator+=(difference_type n) + { + for (int i = 0; i < n; ++i) { - return !(*this == other); + (*this)++; } + return *this; + } - private: - pointer data_; - DataLayout layout_; - int start_; - int end_; - }; +template +typename DataView::Iterator DataView::Iterator::operator+(difference_type n) + { + Iterator tmp(*this); + tmp += n; + return tmp; + } + +template +bool DataView::Iterator::operator==(const Iterator& other) const + { + return (get() == other.get()); + } + +template +bool DataView::Iterator::operator!=(const Iterator& other) const + { + return !(*this == other); + } } // namespace flyft diff --git a/include/flyft/field.h b/include/flyft/field.h index a32e974..b24b841 100644 --- a/include/flyft/field.h +++ b/include/flyft/field.h @@ -1,165 +1,17 @@ #ifndef FLYFT_FIELD_H_ #define FLYFT_FIELD_H_ -#include "flyft/data_layout.h" -#include "flyft/data_view.h" -#include "flyft/tracked_object.h" - -#include -#include +#include "flyft/data_array.h" namespace flyft { -template -class GenericField : public TrackedObject - { - public: - using View = DataView; - using ConstantView = DataView; - using Iterator = typename View::Iterator; - using ConstantIterator = typename ConstantView::Iterator; - - GenericField() = delete; - explicit GenericField(int shape) : GenericField(shape, 0) {} - GenericField(int shape, int buffer_shape) - : data_(nullptr), shape_(0), buffer_shape_(0), layout_(0) - { - reshape(shape, buffer_shape); - } - - // noncopyable / nonmovable - GenericField(const GenericField&) = delete; - GenericField(GenericField&&) = delete; - GenericField& operator=(const GenericField&) = delete; - GenericField& operator=(GenericField&&) = delete; - - ~GenericField() - { - if (data_) - delete[] data_; - } - - T& operator()(int idx) - { - return view()(idx); - } - - const T& operator()(int idx) const - { - return const_view()(idx); - } - - int shape() const - { - return shape_; - } - - int buffer_shape() const - { - return buffer_shape_; - } - - int full_shape() const - { - return shape_ + 2 * buffer_shape_; - } - - View view() - { - token_.stageAndCommit(); - return View(data_, layout_, buffer_shape_, full_shape() - buffer_shape_); - } - - ConstantView const_view() const - { - return ConstantView(static_cast(data_), - layout_, - buffer_shape_, - full_shape() - buffer_shape_); - } - - View full_view() - { - token_.stageAndCommit(); - return View(data_, layout_, 0, full_shape()); - } - - ConstantView const_full_view() const - { - return ConstantView(data_, layout_, 0, full_shape()); - } - - void reshape(int shape, int buffer_shape) - { - if (shape != shape_ || buffer_shape != buffer_shape_) - { - DataLayout layout(shape + 2 * buffer_shape); - - // if data is alloc'd, decide what to do with it - if (data_ != nullptr) - { - if (shape == shape_) - { - // data shape is the same but buffer changed, copy - T* tmp = new T[layout.size()]; - std::fill(tmp, tmp + layout.size(), T(0)); - - // copy the data from old to new - auto view_old = const_view(); - DataView view_new(tmp, layout, buffer_shape, buffer_shape + shape); - std::copy(view_old.begin(), view_old.end(), view_new.begin()); - std::swap(data_, tmp); - delete[] tmp; - } - else if (layout.size() == layout_.size() && layout.size() > 0) - { - // total shape is still the same, just reset the data - std::fill(data_, data_ + layout.size(), T(0)); - } - else - { - // shape is different, wipe it out and realloc - delete[] data_; - data_ = nullptr; - } - } - - // if data is still not alloc'd, make it so - if (data_ == nullptr && layout.size() > 0) - { - data_ = new T[layout.size()]; - std::fill(data_, data_ + layout.size(), T(0)); - } - shape_ = shape; - buffer_shape_ = buffer_shape; - layout_ = layout; - token_.stageAndCommit(); - } - } - - void setBuffer(int buffer_shape) - { - reshape(shape_, buffer_shape); - } - - void requestBuffer(int buffer_shape) - { - if (buffer_shape > buffer_shape_) - { - setBuffer(buffer_shape); - } - } - - private: - T* data_; - int shape_; - int buffer_shape_; - DataLayout layout_; - }; - -using Field = GenericField; -using ComplexField = GenericField>; +//! Real-valued field. +/*! + * A Field represents a real-valued scalar quantity discretized onto a coordinate space. It is a + * multidimensional DataArray and requires another object to define the coordinates. + */ +using Field = DataArray; } // namespace flyft diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h index ab85cc0..10a6879 100644 --- a/include/flyft/mesh.h +++ b/include/flyft/mesh.h @@ -8,96 +8,322 @@ namespace flyft { - +//! Mesh +/*! + * A mesh is a discretized coordinate space. It has one position coordinate and may optionally have + * one, two, or three orientation coordinates. The position coordinate is a symmetry of a standard + * three-dimensional coordinate system (e.g., Cartesian or spherical). + * + * The orientation coordinates are body-fixed Euler angles defined using the Z-X-Z convention for + * the body's principal axes. The bounds of these angles are [0, 2pi), [0, pi], and [0, 2pi], + * respectively. + */ class Mesh { public: + // No default constructor Mesh() = delete; - Mesh(double lower_bound, + + //! Constructor + /*! + * \param num_orientation_dim Number of orientation coordinates. + * \param shape Number of cells along each dimension. + * \param lower_bound Lower bound of the mesh. + * \param upper_bound Upper bound of the mesh. + * \param lower_bc Boundary condition at position coordinate lower bound. + * \param upper_bc Boundary condition at position coordinate upper bound. + */ + Mesh(int num_orientation_dim, + int* shape, + double lower_bound, double upper_bound, - int shape, BoundaryType lower_bc, BoundaryType upper_bc); - std::shared_ptr slice(int start, int end) const; + //! Copy constructor. + Mesh(const Mesh& other); - //! Get position on the mesh, defined as center of bin - double center(int i) const; + //! Copy assignment operator. + Mesh& operator=(const Mesh& other); - //! Lower bound of entire mesh - double lower_bound() const; + //! Move constructor. + Mesh(const Mesh&& other); - //! Get lower bound of bin - double lower_bound(int i) const; + //! Move assignment operator. + Mesh& operator=(const Mesh&& other); - //! Upper bound of entire mesh - double upper_bound() const; + //! Destructor + virtual ~Mesh(); - //! Get upper bound of bin - double upper_bound(int i) const; + //! Slice mesh by position coordinate. + /*! + * \param start First position index to include. + * \param end First position index to exclude. + * \return Slice of original mesh. + * + * A mesh may not be sliced more than once. + */ + std::shared_ptr slice(int start, int end) const; + //! Calculate lower bound of cell. + /*! + * \param coordinate Lower bound of cell. + * \param cell Cell multi-index. + */ + void cell_lower_bound(double* coordinate, const int* cell) const; + + //! Calculate lower bound of cell in a given dimension. + /*! + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns Lower bound of cell in specified dimension. + */ + double cell_lower_bound(int dimension, int index) const; + + //! Calculate lower bound of cell in position space. + /*! + * \param coordinate Lower bound of cell in position space. + * \param cell Cell multi-index. + */ + void cell_position_lower_bound(double* coordinate, const int* cell) const; + + //! Calculate lower bound of cell in orientation space. + /*! + * \param coordinate Lower bound of cell in orientation space. + * \param cell Cell multi-index. + */ + void cell_orientation_lower_bound(double* coordinate, const int* cell) const; + + //! Calculate center of cell. + /*! + * \param coordinate Center of cell. + * \param cell Cell multi-index. + */ + void cell_center(double* coordinate, const int* cell) const; + + //! Calculate center of cell in a given dimension. + /*! + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns Center of cell in specified dimension. + */ + double cell_center(int dimension, int index) const; + + //! Calculate center of cell in position space. + /*! + * \param coordinate Center of cell in position space. + * \param cell Cell multi-index. + */ + void cell_position_center(double* coordinate, const int* cell) const; + + //! Calculate center of cell in orientation space. + /*! + * \param coordinate Center of cell in orientation space. + * \param cell Cell multi-index. + */ + void cell_orientation_center(double* coordinate, const int* cell) const; + + //! Calculate upper bound of cell. + /*! + * \param coordinate Upper bound of cell. + * \param cell Cell multi-index. + */ + void cell_upper_bound(double* coordinate, const int* cell) const; + + //! Calculate upper bound of cell in a given dimension. + /*! + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns upper bound of cell in specified dimension. + */ + double cell_upper_bound(int dimension, int index) const; + + //! Calculate upper bound of cell in position space. + /*! + * \param coordinate Upper bound of cell in position space. + * \param cell Cell multi-index. + */ + void cell_position_upper_bound(double* coordinate, const int* cell) const; + + //! Calculate upper bound of cell in orientation space. + /*! + * \param coordinate Upper bound of cell in orientation space. + * \param cell Cell multi-index. + */ + void cell_orientation_upper_bound(double* coordinate, const int* cell) const; + +#if 0 //! Get surface area of lower edge of bin + /*! + * \param i Multidimensional bin index + * \return Upper bound of the multidimensional bin + */ virtual double area(int i) const = 0; - - //! Get volume of mesh - virtual double volume() const = 0; - - //! Get volume of bin - virtual double volume(int i) const = 0; - - //! Get the bin for a coordinate - int bin(double x) const; - - //! Length of the mesh - double L() const; - - //! Shape of the mesh - int shape() const; - - //! Step size of the mesh - double step() const; - - //! Boundary condition on lower bound of mesh +#endif + + //! Calculate volume of cell. + /*! + * \param cell Cell multi-index. + * \return Volume of cell. + * + * When orientation coordinates are present, the cell volume is a hypervolume over position and + * orientation space. + */ + double cell_volume(const int* cell) const; + + //! Calculate volume of cell in position space. + /*! + * \param cell Cell multi-index. + * \return Volume of cell in position space. + * + * The position coordinates reflect symmetries of the three-dimensional coordinate system, so + * this quantity is always a volume regardless of the number of coordinates. + */ + virtual double cell_position_volume(const int* cell) const = 0; + + //! Calculate volume of cell in orientation space. + /*! + * \param cell Cell multi-index. + * \return Volume of cell in orientation space. + * + * Orientation coordinates are only present as required, so this volume is actually a solid + * angle over valid coordinates. It is zero if there are no orientation coordinates. + */ + double cell_orientation_volume(const int* cell) const; + + //! Calculate the cell that contains a coordinate. + /*! + * \param cell Cell multi-index. + * \param coordinate Coordinate. + */ + void compute_cell(int* cell, const double* coordinate) const; + + //! Convert a length to a number of cells. + /*! + * \param shape Number of cells per dimension. + * \param length Length per dimension. + */ + void length_as_shape(int* shape, const double* length) const; + + //! Convert a shape to a length. + /*! + * \param length Length per dimension. + * \param shape Number of cells per dimension. + */ + void shape_as_length(double* length, const int* shape) const; + + //! Number of position coordinates. + int num_position_coordinates() const; + + //! Number of orientation coordinates. + int num_orientation_coordinates() const; + + //! Total number of coordinates. + int num_coordinates() const; + + //! Number of cells per dimension. + const int* shape() const; + + //! Size of cell per dimension. + const double* step() const; + + //! Boundary condition at position coordinate lower bound. BoundaryType lower_boundary_condition() const; - //! Boundary condition on upper bound of mesh + //! Boundary condition at position coordinate upper bound. BoundaryType upper_boundary_condition() const; - int asShape(double dx) const; - - double asLength(int shape) const; - - double integrateSurface(int idx, double j_lo, double j_hi) const; - double integrateSurface(int idx, const DataView& j) const; - double integrateSurface(int idx, const DataView& j) const; +#if 0 + //! Get the integral of the surface over the mesh + /*! + * \param shape Shape of the mesh + * \return Length of the mesh + */ + double integrateSurface(const std::vector& idx, double j_lo, double j_hi) const; + double integrateSurface(const std::vector& idx, const DataView& j) const; + double integrateSurface(const std::vector& idx, const DataView& j) const; +#endif + + //! Integrate function over cell volume. + /*! + * \param cell Cell multi-index. + * \param f Discretized function to integrate. + * \return Integral of the function over the cell volume. + * + * The integration is performed using a simple midpoint method. + */ + template + double integrate_cell_volume(const int* cell, const T& f) const + { + return cell_volume(cell) * f(cell); + } - double integrateVolume(int idx, double f) const; - double integrateVolume(int idx, const DataView& f) const; - double integrateVolume(int idx, const DataView& f) const; + //! Integrate function over cell orientation volume. + /*! + * \param cell Cell multi-index. + * \param f Discretized function to integrate. + * \return Integral of the function over the cell orientation volume. + * + * The integration is performed using a simple midpoint method. + */ + template + double integrate_cell_orientation_volume(const int* cell, const T& f) const + { + return cell_orientation_volume(cell) * f(cell); + } +#if 0 + //! Interpolate a function on the mesh + /*! + * \param x Coordinate position of the bin + * \param f Function to interpolate + * \return Interpolated function value + * + * The function is interpolated using linear interpolation. + */ template typename std::remove_const::type interpolate(double x, const DataView& f) const; - virtual double gradient(int idx, double f_lo, double f_hi) const = 0; - double gradient(int idx, const DataView& f) const; - double gradient(int idx, const DataView& f) const; - + //! Get the gradient of the mesh + virtual double gradient(const std::vector& idx, double f_lo, double f_hi) const = 0; + double gradient(const std::vector& idx, const DataView& f) const; + double gradient(const std::vector& idx, const DataView& f) const; +#endif + + //! Check if two meshes are equivalent. + /*! + * The meshes are identical if they have the same type, same numbers of coordinates, same + * coordinate space, same shape, and same boundary conditions. + */ bool operator==(const Mesh& other) const; + + //! Check if two meshes are not equivalent. bool operator!=(const Mesh& other) const; + //! Check if boundary conditions are valid. + virtual bool validate_boundary_conditions() const; + protected: - double lower_; - int shape_; //!< Shape of the mesh - BoundaryType lower_bc_; - BoundaryType upper_bc_; - double step_; //!< Spacing between mesh points - int start_; + const int num_position_dim_ = 1; //!< Number of position coordinates + int num_orientation_dim_; //!< Number of orientation coordinates + int num_dim_; //!< Total number of coordinates + int* shape_; //!< Shape of the mesh + int start_; //!< First index in mesh - void validateBoundaryCondition() const; + double* lower_; //!< Lower bounds of coordinates + double* step_; //!< Spacing between mesh points + + BoundaryType lower_bc_; //!< Boundary condition of the lower bound + BoundaryType upper_bc_; //!< Boundary condition of the upper bound virtual std::shared_ptr clone() const = 0; + + private: + //! Allocate per dimension memory. + void allocate(); }; +#if 0 template typename std::remove_const::type Mesh::interpolate(double x, const DataView& f) const { @@ -122,6 +348,7 @@ typename std::remove_const::type Mesh::interpolate(double x, const DataView local() const; std::shared_ptr full() const; + //! Number of position coordinates. + int num_position_coordinates() const; + + //! Number of orientation coordinates. + int num_orientation_coordinates() const; + + //! Total number of coordinates. + int num_coordinates() const; + int getProcessorShape() const; int getProcessorCoordinates() const; int getProcessorCoordinatesByOffset(int offset) const; - int findProcessor(int idx) const; + int findProcessor(const int* cell) const; void sync(std::shared_ptr field); void startSync(std::shared_ptr field); @@ -38,12 +47,11 @@ class ParallelMesh std::shared_ptr gather(std::shared_ptr field, int root) const; private: - std::shared_ptr comm_; - DataLayout layout_; - int coords_; + std::shared_ptr comm_; //!< Communicator std::shared_ptr full_mesh_; std::shared_ptr local_mesh_; + std::vector starts_; std::vector ends_; diff --git a/include/flyft/spherical_mesh.h b/include/flyft/spherical_mesh.h index 9583ce5..f0dd76e 100644 --- a/include/flyft/spherical_mesh.h +++ b/include/flyft/spherical_mesh.h @@ -10,19 +10,13 @@ namespace flyft class SphericalMesh : public Mesh { public: - SphericalMesh(double lower_bound, - double upper_bound, - int shape, - BoundaryType lower_bc, - BoundaryType upper_bc); + // double area(int i) const override; + double cell_position_volume(const int* cell) const override; + // double gradient(int idx, double f_lo, double f_hi) const override; - double area(int i) const override; - double volume() const override; - double volume(int i) const override; - double gradient(int idx, double f_lo, double f_hi) const override; + bool validate_boundary_conditions() const override; protected: - void validateBoundaryCondition() const; std::shared_ptr clone() const override; }; } // namespace flyft diff --git a/include/flyft/state.h b/include/flyft/state.h index 0d61fe1..51346aa 100644 --- a/include/flyft/state.h +++ b/include/flyft/state.h @@ -15,60 +15,124 @@ namespace flyft { +//! Physical state +/*! + * The State defines the density fields that can be operated on. It consists of a mesh that + * discretizes space (and may be parallelized to multiple processors) and one density field on that + * mesh for each type in its specified list. The state also has a time associated with it, which may + * be used for advanced by dynamic calculations. + */ class State : public TrackedObject { public: State() = delete; + + //! Constructor (single type). + /*! + * \param mesh Mesh defining discretized position and orientation coordinates. + * \param type Single type for density field. + */ State(std::shared_ptr mesh, const std::string& type); + + //! Constructor (multiple types). + /*! + * \param mesh Mesh defining discretized position and orientation coordinates. + * \param type List of types for density fields. + */ State(std::shared_ptr mesh, const std::vector& types); + + //! Copy constructor. State(const State& other); + + //! Move constructor. State(State&& other); + + //! Copy assignment operator. State& operator=(const State& other); + + //! Move assignment operator. State& operator=(State&& other); + //! Parallel mesh. std::shared_ptr getMesh(); + + //! Parallel mesh. std::shared_ptr getMesh() const; + + //! Communicator. std::shared_ptr getCommunicator(); + + //! Communicator. std::shared_ptr getCommunicator() const; + //! Number of fields. int getNumFields() const; + + //! List of types. const std::vector& getTypes(); + + //! Get type name associated with type index. const std::string getType(int idx) const; + + //! Get type index associated with type name. int getTypeIndex(const std::string& type) const; + //! Get fields for all types. const TypeMap>& getFields(); + + //! Get field for a single type. std::shared_ptr getField(const std::string& type); + + //! Get field for a single type. std::shared_ptr getField(const std::string& type) const; - void requestFieldBuffer(const std::string& type, int buffer_request); + //! Request a buffer of a particle shape for a field of a single type. + void requestFieldBuffer(const std::string& type, const std::vector& buffer_request); + + //! Gather all fields onto a rank. TypeMap> gatherFields(int rank) const; + + //! Gather a field for a single type onto a rank. std::shared_ptr gatherField(const std::string& type, int rank) const; + //! Synchronize all fields in the state. void syncFields(); + + //! Start synchronizing all fields in the state. + /*! + * This is a nonblocking function, so work can be overlapped with the synchronization. + * Call endSyncFields to finalize. + */ void startSyncFields(); + + //! Finalize sychronization of all fields in the state. void endSyncFields(); - void syncFields(const TypeMap>& fields) const; - void startSyncFields(const TypeMap>& fields) const; - void endSyncFields(const TypeMap>& fields) const; - void endSyncAll() const; + //! Make a set of fields with same types and buffers as the state. void matchFields(TypeMap>& fields) const; + + //! Make a set of fields with same types and buffers as the state, with buffers requested. void matchFields(TypeMap>& fields, - const TypeMap& buffer_requests) const; + const TypeMap>& buffer_requests) const; + //! Get the time. double getTime() const; + + //! Set the time. void setTime(double time); + + //! Advance the time. void advanceTime(double timestep); Token token() override; private: - std::shared_ptr mesh_; - std::vector types_; - TypeMap> fields_; - double time_; + std::shared_ptr mesh_; // types_; //!< Types contained in the state. + TypeMap> fields_; //!< Density fields for the state. + double time_; //!< Time. - Dependencies depends_; + Dependencies depends_; //!< Dependencies on other objects. }; } // namespace flyft diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 84bca5e..f7cd5e8 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,3 +1,4 @@ +find_package(Python REQUIRED COMPONENTS Interpreter Development) find_package(pybind11 CONFIG REQUIRED) add_subdirectory(flyft) if(FLYFT_TESTING) diff --git a/python/flyft/CMakeLists.txt b/python/flyft/CMakeLists.txt index 158ec16..c1fbb2d 100644 --- a/python/flyft/CMakeLists.txt +++ b/python/flyft/CMakeLists.txt @@ -76,11 +76,9 @@ endif() foreach(file IN LISTS FLYFT_PY_SOURCES) add_custom_command( OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${file} - COMMAND ${CMAKE_COMMAND} -E copy - ${CMAKE_CURRENT_LIST_DIR}/${file} - ${CMAKE_CURRENT_BINARY_DIR}/${file} DEPENDS ${file} - POST_BUILD + COMMAND ${CMAKE_COMMAND} + ARGS -E copy ${CMAKE_CURRENT_LIST_DIR}/${file} ${CMAKE_CURRENT_BINARY_DIR}/${file} COMMENT "Copying flyft/${file}" ) endforeach() diff --git a/python/tests/CMakeLists.txt b/python/tests/CMakeLists.txt index d1ab9df..9a6e8bb 100644 --- a/python/tests/CMakeLists.txt +++ b/python/tests/CMakeLists.txt @@ -36,11 +36,9 @@ set(FLYFT_PYTEST_SOURCES foreach(file IN LISTS FLYFT_PYTEST_SOURCES) add_custom_command( OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${file} - COMMAND ${CMAKE_COMMAND} -E copy - ${CMAKE_CURRENT_LIST_DIR}/${file} - ${CMAKE_CURRENT_BINARY_DIR}/${file} DEPENDS ${file} - POST_BUILD + COMMAND ${CMAKE_COMMAND} + ARGS -E copy ${CMAKE_CURRENT_LIST_DIR}/${file} ${CMAKE_CURRENT_BINARY_DIR}/${file} COMMENT "Copying flyft tests/${file}" ) endforeach() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c0fc40f..c21b73c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,41 +1,42 @@ add_library(flyft SHARED - boublik_hard_sphere_functional.cc - brownian_diffusive_flux.cc + # boublik_hard_sphere_functional.cc + # brownian_diffusive_flux.cc cartesian_mesh.cc - composite_external_potential.cc - composite_flux.cc - composite_functional.cc + # composite_external_potential.cc + # composite_flux.cc + # composite_functional.cc communicator.cc - crank_nicolson_integrator.cc + # crank_nicolson_integrator.cc data_layout.cc - explicit_euler_integrator.cc - exponential_wall_potential.cc - external_potential.cc - flux.cc - fourier_transform.cc - functional.cc - grand_potential.cc - grid_interpolator.cc - hard_wall_potential.cc - harmonic_wall_potential.cc - ideal_gas_functional.cc - implicit_euler_integrator.cc - integrator.cc - lennard_jones_93_wall_potential.cc - linear_potential.cc + data_view.cc + # explicit_euler_integrator.cc + # exponential_wall_potential.cc + # external_potential.cc + # flux.cc + # fourier_transform.cc + # functional.cc + # grand_potential.cc + # grid_interpolator.cc + # hard_wall_potential.cc + # harmonic_wall_potential.cc + # ideal_gas_functional.cc + # implicit_euler_integrator.cc + # integrator.cc + # lennard_jones_93_wall_potential.cc + # linear_potential.cc mesh.cc parallel_mesh.cc - picard_iteration.cc - rosenfeld_fmt.cc - rpy_diffusive_flux.cc - solver.cc + # picard_iteration.cc + # rosenfeld_fmt.cc + # rpy_diffusive_flux.cc + # solver.cc spherical_mesh.cc state.cc - three_dimensional_index.cc + # three_dimensional_index.cc tracked_object.cc - virial_expansion.cc - white_bear.cc - white_bear_mark_ii.cc + # virial_expansion.cc + # white_bear.cc + # white_bear_mark_ii.cc ) target_include_directories( diff --git a/src/cartesian_mesh.cc b/src/cartesian_mesh.cc index b1f194f..ce709c5 100644 --- a/src/cartesian_mesh.cc +++ b/src/cartesian_mesh.cc @@ -1,40 +1,68 @@ #include "flyft/cartesian_mesh.h" +#include + namespace flyft { -CartesianMesh::CartesianMesh(double lower_bound, +CartesianMesh::CartesianMesh(int num_orientation_dim, + int* shape, + double lower_bound, double upper_bound, - int shape, BoundaryType lower_bc, BoundaryType upper_bc, double area) - : Mesh(lower_bound, upper_bound, shape, lower_bc, upper_bc), area_(area) + : Mesh(num_orientation_dim, shape, lower_bound, upper_bound, lower_bc, upper_bc), area_(area) { } -std::shared_ptr CartesianMesh::clone() const +CartesianMesh::CartesianMesh(const CartesianMesh& other) : Mesh(other), area_(other.area_) {} + +CartesianMesh& CartesianMesh::operator=(const CartesianMesh& other) { - return std::make_shared(*this); + Mesh::operator=(other); + if (this != &other) + { + area_ = other.area_; + } + return *this; } -double CartesianMesh::area(int /*i*/) const +CartesianMesh::CartesianMesh(const CartesianMesh&& other) + : Mesh(other), area_(std::move(other.area_)) { - return area_; } -double CartesianMesh::volume() const +CartesianMesh& CartesianMesh::operator=(const CartesianMesh&& other) { - return area_ * L(); + Mesh::operator=(other); + area_ = std::move(other.area_); + return *this; } -double CartesianMesh::volume(int /*i*/) const +#if 0 +double CartesianMesh::area(int /*i*/) const { - return area_ * step_; + return area_; } +#endif +double CartesianMesh::cell_position_volume(const int* /*cell*/) const + { + assert(num_position_dim_ == 1); + return area_ * step_[0]; + } + +#if 0 double CartesianMesh::gradient(int /*i*/, double f_lo, double f_hi) const { return (f_hi - f_lo) / step_; } +#endif + +std::shared_ptr CartesianMesh::clone() const + { + return std::make_shared(*this); + } + } // namespace flyft diff --git a/src/data_layout.cc b/src/data_layout.cc index dc58eab..f6ac194 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -1,37 +1,8 @@ #include "flyft/data_layout.h" -#include - namespace flyft { -DataLayout::DataLayout() : shape_(0) {} - -DataLayout::DataLayout(int shape) : shape_(shape) {} - -int DataLayout::operator()(int idx) const - { - return idx; - } - -int DataLayout::shape() const - { - return shape_; - } - -int DataLayout::size() const - { - return shape(); - } - -bool DataLayout::operator==(const DataLayout& other) const - { - return (shape_ == other.shape_); - } - -bool DataLayout::operator!=(const DataLayout& other) const - { - return !(*this == other); - } +template class DataLayout<4>; } // namespace flyft diff --git a/src/data_view.cc b/src/data_view.cc new file mode 100644 index 0000000..9f3d7a1 --- /dev/null +++ b/src/data_view.cc @@ -0,0 +1,9 @@ +#include "flyft/data_view.h" + +namespace flyft + { + +template class DataView; +template class DataView; + + } // namespace flyft diff --git a/src/mesh.cc b/src/mesh.cc index 6e063f3..053f632 100644 --- a/src/mesh.cc +++ b/src/mesh.cc @@ -1,125 +1,334 @@ #include "flyft/mesh.h" +#include #include namespace flyft { - -Mesh::Mesh(double lower_bound, +Mesh::Mesh(int num_orientation_dim, + int* shape, + double lower_bound, double upper_bound, - int shape, BoundaryType lower_bc, BoundaryType upper_bc) - : lower_(lower_bound), shape_(shape), lower_bc_(lower_bc), upper_bc_(upper_bc), start_(0) + : num_orientation_dim_(num_orientation_dim), num_dim_(num_position_dim_ + num_orientation_dim_), + start_(0), lower_bc_(lower_bc), upper_bc_(upper_bc) { - step_ = (upper_bound - lower_bound) / shape_; - validateBoundaryCondition(); + allocate(); + + // copy position information + assert(num_position_dim_ == 1); + shape_[0] = shape[0]; + lower_[0] = lower_bound; + step_[0] = (upper_bound - lower_bound) / shape[0]; + + // copy angle bounds + double angle_lower_bounds[3] {0, 0, 0}; + double angle_upper_bounds[3] {2 * M_PI, M_PI, 2 * M_PI}; + for (int dim = num_position_dim_; dim < num_dim_; ++dim) + { + shape_[dim] = shape_[dim]; + lower_[dim] = angle_lower_bounds[dim - num_position_dim_]; + step_[dim] = (angle_upper_bounds[dim - num_position_dim_] - lower_[dim]) / shape_[dim]; + } } -std::shared_ptr Mesh::slice(int start, int end) const +Mesh::Mesh(const Mesh& other) + : num_orientation_dim_(other.num_orientation_dim_), num_dim_(other.num_dim_), + start_(other.start_), lower_bc_(other.lower_bc_), upper_bc_(other.upper_bc_) { - if (lower_bc_ == BoundaryType::internal || upper_bc_ == BoundaryType::internal) + allocate(); + for (int dim = 0; dim < num_dim_; ++dim) { - throw std::runtime_error("Cannot slice a Mesh more than once."); + shape_[dim] = other.shape_[dim]; + lower_[dim] = other.lower_[dim]; + step_[dim] = other.step_[dim]; } + } + +Mesh& Mesh::operator=(const Mesh& other) + { + if (this != &other) + { + num_orientation_dim_ = other.num_orientation_dim_; + num_dim_ = other.num_dim_; + lower_bc_ = other.lower_bc_; + upper_bc_ = other.upper_bc_; + start_ = other.start_; + + allocate(); + for (int dim = 0; dim < num_dim_; ++dim) + { + shape_[dim] = other.shape_[dim]; + lower_[dim] = other.lower_[dim]; + step_[dim] = other.step_[dim]; + } + } + return *this; + } + +Mesh::Mesh(const Mesh&& other) + : num_orientation_dim_(std::move(other.num_orientation_dim_)), + num_dim_(std::move(other.num_dim_)), shape_(std::move(other.shape_)), + start_(std::move(other.start_)), lower_(std::move(other.lower_)), + step_(std::move(other.step_)), lower_bc_(std::move(other.lower_bc_)), + upper_bc_(std::move(other.upper_bc_)) + + { + } + +Mesh& Mesh::operator=(const Mesh&& other) + { + num_orientation_dim_ = std::move(other.num_orientation_dim_); + num_dim_ = std::move(other.num_dim_); + shape_ = std::move(other.shape_); + start_ = std::move(other.start_); + lower_ = std::move(other.lower_); + step_ = std::move(other.step_); + lower_bc_ = std::move(other.lower_bc_); + upper_bc_ = std::move(other.upper_bc_); + return *this; + } + +Mesh::~Mesh() + { + delete[] shape_; + delete[] lower_; + delete[] step_; + } +std::shared_ptr Mesh::slice(int start, int end) const + { auto m = clone(); + + assert(num_position_dim_ == 1); m->start_ = start; - m->shape_ = end - start; + m->shape_[0] = end - start; if (start > 0) { m->lower_bc_ = BoundaryType::internal; } - if (end < shape_) + if (end < shape_[0]) { m->upper_bc_ = BoundaryType::internal; } + return m; } -double Mesh::center(int i) const +void Mesh::cell_lower_bound(double* coordinate, const int* cell) const { - return lower_ + static_cast(start_ + i + 0.5) * step_; + cell_position_lower_bound(coordinate, cell); + cell_orientation_lower_bound(coordinate + num_position_dim_, cell); } -int Mesh::bin(double x) const +double Mesh::cell_lower_bound(int dimension, int index) const { - return static_cast((x - lower_) / step_) - start_; + if (dimension < num_position_dim_) + { + assert(num_position_dim_ == 1); + return lower_[0] + (start_ + index) * step_[0]; + } + else + { + return lower_[dimension] + index * step_[dimension]; + } } -double Mesh::lower_bound() const +void Mesh::cell_position_lower_bound(double* coordinate, const int* cell) const { - return lower_bound(0); + assert(num_position_dim_ == 1); + coordinate[0] = cell_lower_bound(0, cell[0]); } -double Mesh::lower_bound(int i) const +void Mesh::cell_orientation_lower_bound(double* coordinate, const int* cell) const { - return lower_ + (start_ + i) * step_; + for (int orientation_dim = 0; orientation_dim < num_orientation_dim_; ++orientation_dim) + { + coordinate[orientation_dim] + = cell_lower_bound(num_position_dim_ + orientation_dim, cell[orientation_dim]); + } } -double Mesh::upper_bound() const +void Mesh::cell_center(double* coordinate, const int* cell) const { - return lower_bound(shape_); + cell_position_center(coordinate, cell); + cell_orientation_center(coordinate + num_position_dim_, cell); } -double Mesh::upper_bound(int i) const +double Mesh::cell_center(int dimension, int index) const { - return lower_bound(i + 1); + if (dimension < num_position_dim_) + { + assert(num_position_dim_ == 1); + return lower_[0] + (start_ + index + 0.5) * step_[0]; + } + else + { + return lower_[dimension] + (index + 0.5) * step_[dimension]; + } } -double Mesh::L() const +void Mesh::cell_position_center(double* coordinate, const int* cell) const { - return asLength(shape_); + assert(num_position_dim_ == 1); + coordinate[0] = cell_center(0, cell[0]); } -int Mesh::shape() const +void Mesh::cell_orientation_center(double* coordinate, const int* cell) const { - return shape_; + for (int orientation_dim = 0; orientation_dim < num_orientation_dim_; ++orientation_dim) + { + coordinate[orientation_dim] + = cell_center(num_position_dim_ + orientation_dim, cell[orientation_dim]); + } } -double Mesh::step() const +void Mesh::cell_upper_bound(double* coordinate, const int* cell) const { - return step_; + cell_position_upper_bound(coordinate, cell); + cell_orientation_upper_bound(coordinate + num_position_dim_, cell); } -int Mesh::asShape(double dx) const +double Mesh::cell_upper_bound(int dimension, int index) const { - return static_cast(std::ceil(dx / step_)); + if (dimension < num_position_dim_) + { + assert(num_position_dim_ == 1); + return lower_[0] + (start_ + index + 1) * step_[0]; + } + else + { + return lower_[dimension] + (index + 1) * step_[dimension]; + } } -double Mesh::asLength(int shape) const +void Mesh::cell_position_upper_bound(double* coordinate, const int* cell) const { - return shape * step_; + assert(num_position_dim_ == 1); + coordinate[0] = cell_upper_bound(0, cell[0]); } -double Mesh::integrateSurface(int idx, double j_lo, double j_hi) const +void Mesh::cell_orientation_upper_bound(double* coordinate, const int* cell) const { - return area(idx) * j_lo - area(idx + 1) * j_hi; + for (int orientation_dim = 0; orientation_dim < num_orientation_dim_; ++orientation_dim) + { + coordinate[orientation_dim] + = cell_upper_bound(num_position_dim_ + orientation_dim, cell[orientation_dim]); + } } -double Mesh::integrateSurface(int idx, const DataView& j) const +double Mesh::cell_volume(const int* cell) const { - return integrateSurface(idx, j(idx), j(idx + 1)); + return cell_position_volume(cell) * cell_orientation_volume(cell); } -double Mesh::integrateSurface(int idx, const DataView& j) const +double Mesh::cell_orientation_volume(const int* cell) const { - return integrateSurface(idx, j(idx), j(idx + 1)); + // start with integral being zero in case no orientation coordinates + double V = 0; + + // alpha integral is step size + if (num_orientation_dim_ >= 1) + { + V = step_[num_position_dim_ + 0]; + } + + // beta integral is -cos(beta) + if (num_orientation_dim_ >= 2) + { + const int beta_dim = num_position_dim_ + 1; + const double beta_0 = cell_lower_bound(beta_dim, cell[beta_dim]); + const double beta_1 = cell_upper_bound(beta_dim, cell[beta_dim]); + V *= (std::cos(beta_0) - std::cos(beta_1)); + } + + // gamma integral is step size + if (num_orientation_dim_ == 3) + { + V *= step_[num_position_dim_ + 2]; + } + + return V; + } + +void Mesh::compute_cell(int* cell, const double* coordinate) const + { + for (int dim = 0; dim < num_dim_; ++dim) + { + cell[dim] = static_cast((coordinate[dim] - lower_[dim]) / step_[dim]); + } + + // shift position coordinate + assert(num_position_dim_ == 1); + cell[0] -= start_; + } + +void Mesh::length_as_shape(int* shape, const double* length) const + { + for (int dim = 0; dim < num_dim_; ++dim) + { + shape[dim] = static_cast(std::ceil(length[dim] / step_[dim])); + } + } + +void Mesh::shape_as_length(double* length, const int* shape) const + { + for (int dim = 0; dim < num_dim_; ++dim) + { + length[dim] = shape[dim] * step_[dim]; + } + } + +int Mesh::num_position_coordinates() const + { + return num_position_dim_; + } + +int Mesh::num_orientation_coordinates() const + { + return num_orientation_dim_; + } + +int Mesh::num_coordinates() const + { + return num_dim_; + } + +const int* Mesh::shape() const + { + return shape_; + } + +const double* Mesh::step() const + { + return step_; + } + +BoundaryType Mesh::lower_boundary_condition() const + { + return lower_bc_; + } + +BoundaryType Mesh::upper_boundary_condition() const + { + return upper_bc_; } -double Mesh::integrateVolume(int idx, double f) const +#if 0 +double Mesh::integrateSurface(const std::vector& idx, double j_lo, double j_hi) const { - return volume(idx) * f; + return area(idx) * j_lo - area(idx + 1) * j_hi; } -double Mesh::integrateVolume(int idx, const DataView& f) const +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const { - return integrateVolume(idx, f(idx)); + return integrateSurface(idx, j(idx), j(idx + 1)); } -double Mesh::integrateVolume(int idx, const DataView& f) const +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const { - return integrateVolume(idx, f(idx)); + return integrateSurface(const std::vector& idx, j(idx), j(idx + 1)); } double Mesh::gradient(int idx, const DataView& f) const @@ -147,12 +356,18 @@ double Mesh::gradient(int idx, const DataView& f) const return gradient(idx, f(idx - 1), f(idx)); } } +#endif bool Mesh::operator==(const Mesh& other) const { - return (typeid(*this) == typeid(other) && lower_ == other.lower_ && shape_ == other.shape_ - && lower_bc_ == other.lower_bc_ && upper_bc_ == other.upper_bc_ - && start_ == other.start_); + bool is_same = (typeid(*this) == typeid(other) && num_position_dim_ == other.num_position_dim_ + && num_orientation_dim_ == other.num_orientation_dim_ && start_ == other.start_ + && lower_bc_ == other.lower_bc_ && upper_bc_ == other.upper_bc_); + for (int dim = 0; dim < num_dim_ && is_same; ++dim) + { + is_same |= (shape_[dim] == other.shape_[dim] && lower_[dim] == other.lower_[dim]); + } + return is_same; } bool Mesh::operator!=(const Mesh& other) const @@ -160,24 +375,26 @@ bool Mesh::operator!=(const Mesh& other) const return !(*this == other); } -BoundaryType Mesh::lower_boundary_condition() const +bool Mesh::validate_boundary_conditions() const { - return lower_bc_; + bool invalid + = ((upper_bc_ == BoundaryType::periodic + && !(lower_bc_ == BoundaryType::periodic || lower_bc_ == BoundaryType::internal)) + || (lower_bc_ == BoundaryType::periodic + && !(upper_bc_ == BoundaryType::periodic || upper_bc_ == BoundaryType::internal))); + return !invalid; } -BoundaryType Mesh::upper_boundary_condition() const +void Mesh::allocate() { - return upper_bc_; - } + delete[] shape_; + shape_ = new int[num_dim_]; -void Mesh::validateBoundaryCondition() const - { - if ((upper_bc_ == BoundaryType::periodic - && !(lower_bc_ == BoundaryType::periodic || lower_bc_ == BoundaryType::internal)) - || (lower_bc_ == BoundaryType::periodic - && !(upper_bc_ == BoundaryType::periodic || upper_bc_ == BoundaryType::internal))) - { - throw std::invalid_argument("Both boundaries must be periodic if one is."); - } + delete[] lower_; + lower_ = new double[num_dim_]; + + delete[] step_; + step_ = new double[num_dim_]; } + } // namespace flyft diff --git a/src/parallel_mesh.cc b/src/parallel_mesh.cc index d8f98e1..a2152d3 100644 --- a/src/parallel_mesh.cc +++ b/src/parallel_mesh.cc @@ -1,46 +1,48 @@ #include "flyft/parallel_mesh.h" +#include #include namespace flyft { ParallelMesh::ParallelMesh(std::shared_ptr mesh, std::shared_ptr comm) + : comm_(comm), full_mesh_(mesh) { - // set the *full* mesh passed to the communicator - full_mesh_ = mesh; - - // set the communicator (coordinates are assumed to be the same as rank for now) - comm_ = comm; - layout_ = DataLayout(comm_->size()); - coords_ = comm_->rank(); - - // determine the mesh sites that each processor owns - const int floor_shape = full_mesh_->shape() / layout_.shape(); - const int leftover = full_mesh_->shape() - layout_.shape() * floor_shape; - starts_ = std::vector(layout_.size(), 0); - ends_ = std::vector(layout_.size(), 0); - for (int idx = 0; idx < layout_.shape(); ++idx) + // validate and set the *full* mesh passed to the communicator + assert(full_mesh_->num_position_coordinates() == 1); + if (!full_mesh_->validate_boundary_conditions()) + { + throw std::runtime_error("Invalid boundary conditions for mesh"); + } + + // determine the mesh cells that each processor owns + const auto num_procs = comm_->size(); + const int floor_shape = full_mesh_->shape()[0] / num_procs; + const int leftover = full_mesh_->shape()[0] - num_procs * floor_shape; + starts_ = std::vector(num_procs, 0); + ends_ = std::vector(num_procs, 0); + for (int idx = 0; idx < num_procs; ++idx) { - const int coord_idx = layout_(idx); if (idx < leftover) { - starts_[coord_idx] = idx * floor_shape + idx; - ends_[coord_idx] = starts_[coord_idx] + (floor_shape + 1); + starts_[idx] = idx * floor_shape + idx; + ends_[idx] = starts_[idx] + (floor_shape + 1); } else { - starts_[coord_idx] = idx * floor_shape + leftover; - ends_[coord_idx] = starts_[coord_idx] + floor_shape; + starts_[idx] = idx * floor_shape + leftover; + ends_[idx] = starts_[idx] + floor_shape; } } // size the local mesh based on the sites covered - const int start = starts_[layout_(coords_)]; - const int end = ends_[layout_(coords_)]; + const auto rank = comm_->rank(); + const int start = starts_[rank]; + const int end = ends_[rank]; if (start == end) { - // ERROR: cannot have empty processor + throw std::runtime_error("Empty processor sized"); } local_mesh_ = full_mesh_->slice(start, end); @@ -68,42 +70,60 @@ std::shared_ptr ParallelMesh::full() const return full_mesh_; } +int ParallelMesh::num_position_coordinates() const + { + return full_mesh_->num_position_coordinates(); + } + +int ParallelMesh::num_orientation_coordinates() const + { + return full_mesh_->num_orientation_coordinates(); + } + +int ParallelMesh::num_coordinates() const + { + return full_mesh_->num_coordinates(); + } + int ParallelMesh::getProcessorShape() const { - return layout_.shape(); + return comm_->size(); } int ParallelMesh::getProcessorCoordinates() const { - return coords_; + return comm_->rank(); } int ParallelMesh::getProcessorCoordinatesByOffset(int offset) const { - int proc = coords_ + offset; + int proc = comm_->rank() + offset; + const auto num_ranks = comm_->size(); while (proc < 0) { - proc += layout_.shape(); + proc += num_ranks; } - while (proc >= layout_.shape()) + while (proc >= num_ranks) { - proc -= layout_.shape(); + proc -= num_ranks; } + assert(proc >= 0 && proc < num_ranks); return proc; } -int ParallelMesh::findProcessor(int idx) const +int ParallelMesh::findProcessor(const int* cell) const { - if (idx < 0 || idx >= full_mesh_->shape()) + if (cell[0] < 0 || cell[0] >= full_mesh_->shape()[0]) { - // ERROR: processor out of range + throw std::runtime_error("Cell out of range"); } // TODO: replace with binary search since vectors are sorted + const auto num_ranks = comm_->size(); int proc = -1; - for (int i = 0; proc < 0 && i < layout_.shape(); ++i) + for (int i = 0; proc < 0 && i < num_ranks; ++i) { - if (idx >= starts_[layout_(i)] && idx < ends_[layout_(i)]) + if (cell[0] >= starts_[i] && cell[0] < ends_[i]) { proc = i; } @@ -135,67 +155,133 @@ void ParallelMesh::startSync(std::shared_ptr field) } #endif - // check field shape - const int shape = field->shape(); - const int buffer_shape = field->buffer_shape(); + // check field shape in decomposed dimension + const int shape = field->shape()[0]; + const int buffer_shape = field->buffer_shape()[0]; if (buffer_shape > shape) { - // ERROR: overdecomposed (only nearest-neighbor comms supported) - throw std::runtime_error("Mesh overdecomposed"); + throw std::runtime_error( + "Mesh overdecomposed, only nearest neighbor communication supported."); } else if (buffer_shape == 0) { // nothing to do, no buffer needed return; } + // sync field auto f = field->view(); const auto lower_bc = local_mesh_->lower_boundary_condition(); const auto upper_bc = local_mesh_->upper_boundary_condition(); #ifdef FLYFT_MPI - const int left = layout_(getProcessorCoordinatesByOffset(-1)); - const int right = layout_(getProcessorCoordinatesByOffset(1)); + const int left = getProcessorCoordinatesByOffset(-1); + const int right = getProcessorCoordinatesByOffset(1); std::vector requests; requests.reserve(4); + // stride for sending data of decomposed coordinate + size_t decomp_stride = 1; + if (comm_->size() > 1 && (lower_bc == BoundaryType::periodic || lower_bc == BoundaryType::internal)) { + // calculate stride per position coordinate + const auto orientation_shape = full_mesh_->shape() + full_mesh_->num_position_coordinates(); + for (int dim = 0; dim < full_mesh_->num_orientation_coordinates(); ++dim) + { + decomp_stride *= orientation_shape[dim]; + } + + // get iterators (pointers) to portions of data we will send / receive from + std::vector multi_index(f.num_dimensions(), 0); + multi_index[0] = -buffer_shape; + const auto recv_it = f.make_iterator(multi_index.data()); + multi_index[0] = 0; + const auto send_it = f.make_iterator(multi_index.data()); + // receive left buffer from left (tag 0), right buffer from right (tag 1) MPI_Comm comm = comm_->get(); const auto end = requests.size(); requests.resize(end + 2); - MPI_Irecv(&f(-buffer_shape), buffer_shape, MPI_DOUBLE, left, 0, comm, &requests[end]); - MPI_Isend(&f(0), buffer_shape, MPI_DOUBLE, left, 1, comm, &requests[end + 1]); + MPI_Irecv(recv_it.get(), + buffer_shape * decomp_stride, + MPI_DOUBLE, + left, + 0, + comm, + &requests[end]); + MPI_Isend(send_it.get(), + buffer_shape * decomp_stride, + MPI_DOUBLE, + left, + 1, + comm, + &requests[end + 1]); } else #endif { - for (int idx = 0; idx < buffer_shape; ++idx) + std::vector multi_index(f.num_dimensions(), 0); + if (lower_bc == BoundaryType::zero) { - double value; - if (lower_bc == BoundaryType::zero) - { - value = 0; - } - else if (lower_bc == BoundaryType::repeat) - { - value = f(0); - } - else if (lower_bc == BoundaryType::reflect) + multi_index[0] = -buffer_shape; + auto recv_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = 0; + const auto recv_it_end = f.make_iterator(multi_index.data()); + + while (recv_it_start != recv_it_end) { - value = f(1 + idx); + *recv_it_start = 0; + ++recv_it_start; } - else if (lower_bc == BoundaryType::periodic || lower_bc == BoundaryType::internal) + } + else if (lower_bc == BoundaryType::repeat || lower_bc == BoundaryType::reflect) + { + for (int idx = 0; idx < buffer_shape; ++idx) { - value = f(shape - 1 - idx); + multi_index[0] = -1 - idx; + auto recv_it_start = f.make_iterator(multi_index.data()); + + ++multi_index[0]; + const auto recv_it_end = f.make_iterator(multi_index.data()); + + if (lower_bc == BoundaryType::repeat) + { + multi_index[0] = 0; + } + else // reflect + { + multi_index[0] = 1 + idx; + } + auto send_it_start = f.make_iterator(multi_index.data()); + + while (recv_it_start != recv_it_end) + { + *recv_it_start = *send_it_start; + ++recv_it_start; + ++send_it_start; + } } - else + } + else // periodic / internal + { + multi_index[0] = -buffer_shape; + auto recv_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = shape - buffer_shape; + auto send_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = shape; + const auto send_it_end = f.make_iterator(multi_index.data()); + + while (send_it_start != send_it_end) { - throw std::runtime_error("Unknown boundary condition"); + *recv_it_start = *send_it_start; + ++recv_it_start; + ++send_it_start; } - f(-1 - idx) = value; } } @@ -203,13 +289,26 @@ void ParallelMesh::startSync(std::shared_ptr field) if (comm_->size() > 1 && (upper_bc == BoundaryType::periodic || upper_bc == BoundaryType::internal)) { + // get iterators (pointers) to portions of data we will send / receive from + std::vector multi_index(f.num_dimensions(), 0); + multi_index[0] = shape; + const auto recv_it = f.make_iterator(multi_index.data()); + multi_index[0] = shape - buffer_shape; + const auto send_it = f.make_iterator(multi_index.data()); + // send left edge to left (tag 1), right edge to right (tag 0) MPI_Comm comm = comm_->get(); const auto end = requests.size(); requests.resize(end + 2); - MPI_Irecv(&f(shape), buffer_shape, MPI_DOUBLE, right, 1, comm, &requests[end]); - MPI_Isend(&f(shape - buffer_shape), - buffer_shape, + MPI_Irecv(recv_it.get(), + buffer_shape * decomp_stride, + MPI_DOUBLE, + right, + 1, + comm, + &requests[end]); + MPI_Isend(send_it.get(), + buffer_shape * decomp_stride, MPI_DOUBLE, right, 0, @@ -219,30 +318,66 @@ void ParallelMesh::startSync(std::shared_ptr field) else #endif { - for (int idx = 0; idx < buffer_shape; ++idx) + std::vector multi_index(f.num_dimensions(), 0); + if (lower_bc == BoundaryType::zero) { - double value; - if (upper_bc == BoundaryType::zero) - { - value = 0; - } - else if (upper_bc == BoundaryType::repeat) - { - value = f(shape - 1); - } - else if (upper_bc == BoundaryType::reflect) + multi_index[0] = shape; + auto recv_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = shape + buffer_shape; + const auto recv_it_end = f.make_iterator(multi_index.data()); + + while (recv_it_start != recv_it_end) { - value = f(shape - 1 - idx); + *recv_it_start = 0; + ++recv_it_start; } - else if (upper_bc == BoundaryType::periodic || upper_bc == BoundaryType::internal) + } + else if (lower_bc == BoundaryType::repeat || lower_bc == BoundaryType::reflect) + { + for (int idx = 0; idx < buffer_shape; ++idx) { - value = f(idx); + multi_index[0] = shape + idx; + auto recv_it_start = f.make_iterator(multi_index.data()); + + ++multi_index[0]; + const auto recv_it_end = f.make_iterator(multi_index.data()); + + if (lower_bc == BoundaryType::repeat) + { + multi_index[0] = shape - 1; + } + else // reflect + { + multi_index[0] = shape - 1 - idx; + } + auto send_it_start = f.make_iterator(multi_index.data()); + + while (recv_it_start != recv_it_end) + { + *recv_it_start = *send_it_start; + ++recv_it_start; + ++send_it_start; + } } - else + } + else // periodic + { + multi_index[0] = shape; + auto recv_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = -buffer_shape; + auto send_it_start = f.make_iterator(multi_index.data()); + + multi_index[0] = 0; + const auto send_it_end = f.make_iterator(multi_index.data()); + + while (send_it_start != send_it_end) { - throw std::runtime_error("Unknown boundary condition"); + *recv_it_start = *send_it_start; + ++send_it_start; + ++recv_it_start; } - f(shape + idx) = value; } } @@ -297,14 +432,19 @@ std::shared_ptr ParallelMesh::gather(std::shared_ptr field, int /* #ifdef FLYFT_MPI if (comm_->size() > 1) { - // determine number of elements sent by each rank - std::vector counts(comm_->size()); - for (int idx = 0; idx < layout_.shape(); ++idx) + // calculate stride per position coordinate + size_t decomp_stride = 1; + const auto orientation_shape = full_mesh_->shape() + full_mesh_->num_position_coordinates(); + for (int dim = 0; dim < full_mesh_->num_orientation_coordinates(); ++dim) { - const auto coord_idx = layout_(idx); - counts[coord_idx] = ends_[coord_idx] - starts_[coord_idx]; + decomp_stride *= orientation_shape[dim]; } + // determine number of elements sent by each rank + const auto num_procs = comm_->size(); + std::vector counts; + std::vector displs; + // send buffer is valid on all ranks const auto f = field->const_view(); @@ -313,18 +453,28 @@ std::shared_ptr ParallelMesh::gather(std::shared_ptr field, int /* void* recv(nullptr); if (comm_->rank() == root) { - new_field = std::make_shared(full_mesh_->shape()); - auto new_f = new_field->view(); - recv = static_cast(&new_f(0)); + new_field = std::make_shared(full_mesh_->num_coordinates(), full_mesh_->shape()); + + recv = static_cast(new_field->view().begin().get()); + counts.resize(num_procs); + displs.resize(num_procs); + for (int idx = 0; idx < num_procs; ++idx) + { + counts[idx] = (ends_[idx] - starts_[idx]) * decomp_stride; + if (idx > 0) + { + displs[idx] = displs[idx - 1] + counts[idx]; + } + } } // gather to the root rank - MPI_Gatherv(&f(0), + MPI_Gatherv(f.begin().get(), f.size(), MPI_DOUBLE, recv, - &counts[0], - &starts_[0], + counts.data(), + displs.data(), MPI_DOUBLE, root, comm_->get()); diff --git a/src/spherical_mesh.cc b/src/spherical_mesh.cc index 698f064..95510fe 100644 --- a/src/spherical_mesh.cc +++ b/src/spherical_mesh.cc @@ -5,57 +5,41 @@ namespace flyft { -SphericalMesh::SphericalMesh(double lower_bound, - double upper_bound, - int shape, - BoundaryType lower_bc, - BoundaryType upper_bc) - : Mesh(lower_bound, upper_bound, shape, lower_bc, upper_bc) - { - validateBoundaryCondition(); - } std::shared_ptr SphericalMesh::clone() const { return std::make_shared(*this); } +#if 0 double SphericalMesh::area(int i) const { const double r = lower_bound(i); return 4. * M_PI * r * r; } +#endif -double SphericalMesh::volume() const +double SphericalMesh::cell_position_volume(const int* cell) const { - const double rlo = lower_bound(); - const double rhi = upper_bound(); - return (4. * M_PI / 3.) * (rhi * rhi * rhi - rlo * rlo * rlo); - } - -double SphericalMesh::volume(int i) const - { - const double r_out = upper_bound(i); - const double r_in = lower_bound(i); + const double r_in = cell_lower_bound(0, cell[0]); + const double r_out = cell_upper_bound(0, cell[0]); return (4. * M_PI / 3.) * (r_out * r_out * r_out - r_in * r_in * r_in); } +#if 0 double SphericalMesh::gradient(int /*i*/, double f_lo, double f_hi) const { return (f_hi - f_lo) / (step_); } +#endif -void SphericalMesh::validateBoundaryCondition() const +bool SphericalMesh::validate_boundary_conditions() const { - Mesh::validateBoundaryCondition(); - - if (lower_bc_ == BoundaryType::periodic || upper_bc_ == BoundaryType::periodic) - { - throw std::invalid_argument("Periodic boundary conditions invalid in spherical geometry"); - } - else if (upper_bc_ == BoundaryType::reflect) - { - throw std::invalid_argument("Reflect boundary condition invalid for upper boundary"); - } + bool invalid = !Mesh::validate_boundary_conditions(); + + invalid |= (lower_bc_ == BoundaryType::periodic || upper_bc_ == BoundaryType::periodic); + invalid |= (upper_bc_ == BoundaryType::reflect); + + return !invalid; } } // namespace flyft diff --git a/src/state.cc b/src/state.cc index 60565cd..161a908 100644 --- a/src/state.cc +++ b/src/state.cc @@ -2,6 +2,7 @@ #include #include +#include namespace flyft { @@ -16,7 +17,7 @@ State::State(std::shared_ptr mesh, const std::vector& { for (const auto& t : types_) { - fields_[t] = std::make_shared(mesh_->local()->shape()); + fields_[t] = std::make_shared(mesh_->num_coordinates(), mesh_->local()->shape()); depends_.add(fields_[t].get()); } } @@ -27,7 +28,9 @@ State::State(const State& other) for (const auto& t : types_) { auto other_field = other.fields_(t); - fields_[t] = std::make_shared(other_field->shape(), other_field->buffer_shape()); + fields_[t] = std::make_shared(other_field->num_dimensions(), + other_field->shape(), + other_field->buffer_shape()); std::copy(other_field->const_full_view().begin(), other_field->const_full_view().end(), fields_[t]->full_view().begin()); @@ -51,10 +54,13 @@ State& State::operator=(const State& other) time_ = other.time_; // match field types and buffer shapes - TypeMap buffers; + TypeMap> buffers; for (const auto& t : types_) { - buffers[t] = other.fields_(t)->buffer_shape(); + const int num_coords = mesh_->full()->num_coordinates(); + const auto other_buffer = other.fields_(t)->buffer_shape(); + buffers[t].resize(num_coords); + std::copy(other_buffer, other_buffer + num_coords, buffers[t].begin()); } other.matchFields(fields_, buffers); @@ -122,7 +128,7 @@ int State::getTypeIndex(const std::string& type) const const auto it = std::find(types_.begin(), types_.end(), type); if (it == types_.end()) { - // error: type not found + throw std::invalid_argument("Type not found"); } return (it - types_.begin()); } @@ -142,9 +148,9 @@ std::shared_ptr State::getField(const std::string& type) const return fields_(type); } -void State::requestFieldBuffer(const std::string& type, int buffer_request) +void State::requestFieldBuffer(const std::string& type, const std::vector& buffer_request) { - fields_(type)->requestBuffer(buffer_request); + fields_(type)->requestBuffer(buffer_request.data()); } TypeMap> State::gatherFields(int rank) const @@ -164,67 +170,33 @@ std::shared_ptr State::gatherField(const std::string& type, int rank) con void State::syncFields() { - syncFields(fields_); + startSyncFields(); + endSyncFields(); } void State::startSyncFields() { - startSyncFields(fields_); - } - -void State::endSyncFields() - { - endSyncFields(fields_); - } - -void State::syncFields(const TypeMap>& fields) const - { - startSyncFields(fields); - endSyncFields(fields); - } - -void State::startSyncFields(const TypeMap>& fields) const - { - for (auto it = fields.cbegin(); it != fields.cend(); ++it) + for (auto it = fields_.cbegin(); it != fields_.cend(); ++it) { - if (std::find(types_.begin(), types_.end(), it->first) == types_.end()) - { - // ERROR: types do not match - } - else - { - mesh_->startSync(it->second); - } + mesh_->startSync(it->second); } } -void State::endSyncFields(const TypeMap>& fields) const +void State::endSyncFields() { - for (auto it = fields.cbegin(); it != fields.cend(); ++it) + for (auto it = fields_.cbegin(); it != fields_.cend(); ++it) { - if (std::find(types_.begin(), types_.end(), it->first) == types_.end()) - { - // ERROR: types do not match - } - else - { - mesh_->endSync(it->second); - } + mesh_->endSync(it->second); } } -void State::endSyncAll() const - { - mesh_->endSyncAll(); - } - void State::matchFields(TypeMap>& fields) const { - matchFields(fields, TypeMap()); + matchFields(fields, TypeMap>()); } void State::matchFields(TypeMap>& fields, - const TypeMap& buffer_requests) const + const TypeMap>& buffer_requests) const { // purge stored types that are not in the state for (auto it = fields.cbegin(); it != fields.cend(); /* no increment here */) @@ -247,7 +219,8 @@ void State::matchFields(TypeMap>& fields, bool has_type = (fields.find(t) != fields.end()); // determine buffer request, preserving buffer if type already exists but no request is made - int buffer_request; + const auto num_dimensions = mesh_->full()->num_coordinates(); + std::vector buffer_request(num_dimensions, 0); auto it = buffer_requests.find(t); if (it != buffer_requests.cend()) { @@ -255,21 +228,20 @@ void State::matchFields(TypeMap>& fields, } else if (has_type) { - buffer_request = fields[t]->buffer_shape(); - } - else - { - buffer_request = 0; + const auto buffer_shape = fields[t]->buffer_shape(); + std::copy(buffer_shape, buffer_shape + num_dimensions, buffer_request.begin()); } // make sure field exists and has the right shape if (!has_type) { - fields[t] = std::make_shared(mesh_->local()->shape(), buffer_request); + fields[t] = std::make_shared(num_dimensions, + mesh_->local()->shape(), + buffer_request.data()); } else { - fields[t]->reshape(mesh_->local()->shape(), buffer_request); + fields[t]->reshape(num_dimensions, mesh_->local()->shape(), buffer_request.data()); } } }