From c2dea16f5ca19ea08bb570f8c759b55055cd6104 Mon Sep 17 00:00:00 2001 From: Mayukh Kundu <79770572+mayukh33@users.noreply.github.com> Date: Wed, 19 Mar 2025 10:58:42 -0500 Subject: [PATCH 01/13] Enable two-dimensional array in the existing data structure (#51) * Add vectors to Data_layout class * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add strided indexing for n-dimensional to 1D array mapping * Add size calculation and updates for doxygen documentation * Fix doxygen comment * Minor fix to data layout comment * Update doxygen documentation * Minor fix to Datalayout class * Apply changes to make start, end, shape and size compatible with multi-dimensional array * Minor fix to data_layout documentation * Fix operator in DataView class * Minor fix * Different implementation of size method in data_layout --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- include/flyft/data_layout.h | 37 +++++++- include/flyft/data_view.h | 90 ++++++++++++++++--- python/tests/conftest.py | 3 +- python/tests/test_boublik_hard_sphere.py | 3 +- python/tests/test_brownian_diffusive_flux.py | 3 +- .../test_composite_external_potential.py | 3 +- python/tests/test_composite_flux.py | 3 +- python/tests/test_composite_functional.py | 3 +- .../tests/test_crank_nicolson_integrator.py | 3 +- .../tests/test_explicit_euler_integrator.py | 3 +- .../tests/test_exponential_wall_potential.py | 3 +- python/tests/test_external_field.py | 3 +- python/tests/test_field.py | 3 +- python/tests/test_grand_potential.py | 3 +- python/tests/test_hard_wall_potential.py | 3 +- python/tests/test_harmonic_wall_potential.py | 3 +- .../tests/test_implicit_euler_integrator.py | 3 +- .../test_lennard_jones_93_wall_potential.py | 3 +- python/tests/test_mirror.py | 3 +- python/tests/test_parameter.py | 3 +- python/tests/test_picard_iteration.py | 3 +- python/tests/test_rpy_diffusive_flux.py | 3 +- python/tests/test_state.py | 3 +- python/tests/test_virial_expansion.py | 3 +- python/tests/test_white_bear.py | 3 +- python/tests/test_white_bear_mark_ii.py | 3 +- python/validate/binary_mixture.py | 3 +- python/validate/brownian_diffusion.py | 3 +- python/validate/hard_sphere_evaporation.py | 3 +- python/validate/hard_sphere_slit.py | 3 +- python/validate/stratification.py | 3 +- src/data_layout.cc | 23 +++-- 32 files changed, 155 insertions(+), 82 deletions(-) diff --git a/include/flyft/data_layout.h b/include/flyft/data_layout.h index ca04fb9..d32077d 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -2,25 +2,54 @@ #define FLYFT_DATA_LAYOUT_H_ #include +#include #include +#include namespace flyft { +//! Multidimensional array layout +/*! + * DataLayout maps an N-dimensional index to a one-dimensional + * index using row-major ordering. Negative indexes are supported + * and interpreted as relative to the zero-index. + */ class DataLayout { public: DataLayout(); - explicit DataLayout(int shape_); + //! Constructor + /*! + * \param shape Shape of the index array. + */ + explicit DataLayout(const std::vector& shape); - int operator()(int idx) const; - int shape() const; + //! Map a multidimensional index to a one-dimensional index + /*! + * \param idx Multidimensional index. + * \return One-dimensional index. + * + * Row-major ordering is used. + */ + int operator()(const std::vector& idx) const; + + //! Shape of the layout in each dimension + std::vector shape() const; + + //! Total number of elements in the layout int size() const; + //! Test if two layouts are equal + /*! + * \return True if the layouts have the same shape. + */ bool operator==(const DataLayout& other) const; + + //! Test if two layouts are not equal bool operator!=(const DataLayout& other) const; private: - int shape_; + std::vector shape_; //!< Multidimensional shape of layout }; } // namespace flyft diff --git a/include/flyft/data_view.h b/include/flyft/data_view.h index 6da3f61..4449d20 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -3,16 +3,23 @@ #include "flyft/data_layout.h" +#include +#include + namespace flyft { template +//! Multidimensional array view +/*! + * DataView provides ability to view a certain section of the multidimensional array + */ 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; // Extracts datatype + using pointer = value_type*; // Pointer to the datatype + using reference = value_type&; // Reference to the datatype class Iterator { @@ -23,6 +30,7 @@ class DataView using pointer = DataView::pointer; using reference = DataView::reference; + // Overloading the constructor of the iterator class Iterator() : Iterator(DataView()) {} explicit Iterator(const DataView& view) : Iterator(view, 0) {} @@ -85,61 +93,117 @@ class DataView int current_; }; + //! Empty constructor + /*! + * The default is a null pointer to a zero-size layout. + */ DataView() : DataView(nullptr, DataLayout()) {} + //! Constructor + /*! + * \param data_ Pointer to the data. + * \param layout_ Layout of the data. + */ + DataView(pointer data, const DataLayout& layout) : DataView(data, layout, 0, layout.shape()) {} - DataView(pointer data, const DataLayout& layout, int start, int end) + //! Constructor + /*! + * \param data_ Pointer to the data. + * \param layout_ Layout of the data. + * \param start_ Start of the array. + * \param end_ End of the array. + * + * Initializes the view of the multi-dimensional array. + */ + DataView(pointer data, const DataLayout& layout, std::vector start, std::vector end) : data_(data), layout_(layout), start_(start), end_(end) { } - reference operator()(int idx) const + //! Reference to the referenced element of the data. + /*! + * \param idx Array of indices for the multiple dimensions + * \return reference to the referenced data + */ + reference operator()(const std::vector& idx) const { - return data_[layout_(start_ + idx)]; + std::vector temp; + std::transform(idx.begin(), idx.end(), start_.begin(), temp.begin(), std::plus()); + return data_[layout_(temp)]; } - int shape() const + //! Shape of the multi-dimensional array + std::vector shape() const { - return end_ - start_; + if (start_.size() != end_.size()) + { + throw std::invalid_argument( + "Arrays must have the same size for element-wise subtraction."); + } + + std::vector& result(start_.size()); + for (int i = 0; i < start_.size(); ++i) + { + result[i] = start_[i] - end_[i]; + } + return result; } + //! Number of elements in the multi-dimensional array int size() const { - return shape(); + return std::accumulate(shape().begin(), shape().end(), 1, std::multiplies()); } + //! Test if pointer to the data is not a null pointer + /*! + * \return DataView is not null if it is a view of valid data. + */ explicit operator bool() const { return (data_ != nullptr); } + //! Start point of the data + /*! + * \return data at index position 0. + */ Iterator begin() const { return Iterator(*this, 0); } + //! End point of the data + /*! + * \return data at end of the array. + */ Iterator end() const { return Iterator(*this, shape()); } + //! Test if two views are equal + /*! + * \return True if the view have the same pointer, layout, start and end points. + */ bool operator==(const DataView& other) const { return (data_ == other.data_ && layout_ == other.layout_ && start_ == other.start_ && end_ == other.end_); } + //! Test if two views are not equal bool operator!=(const DataView& other) const { return !(*this == other); } private: - pointer data_; - DataLayout layout_; - int start_; - int end_; + pointer data_; //!< Pointer to the array + DataLayout layout_; //!< Multidimensional array layout + std::vector start_; //!< Start of the array + std::vector end_; //!< End of the array }; } // namespace flyft diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 04d65ed..4adb701 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -1,8 +1,7 @@ +import flyft import pytest from pytest_lazy_fixtures import lf as lazy_fixture -import flyft - @pytest.fixture def ig(): diff --git a/python/tests/test_boublik_hard_sphere.py b/python/tests/test_boublik_hard_sphere.py index 0b7c3cb..98c4463 100644 --- a/python/tests/test_boublik_hard_sphere.py +++ b/python/tests/test_boublik_hard_sphere.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - def fex_cs(eta, d): """Carnahan-Starling free-energy density of hard spheres""" diff --git a/python/tests/test_brownian_diffusive_flux.py b/python/tests/test_brownian_diffusive_flux.py index f7099c1..e55ad57 100644 --- a/python/tests/test_brownian_diffusive_flux.py +++ b/python/tests/test_brownian_diffusive_flux.py @@ -1,9 +1,8 @@ +import flyft import numpy as np import pytest from pytest_lazy_fixtures import lf as lazy_fixture -import flyft - @pytest.fixture def cartesian_mesh_grand(): diff --git a/python/tests/test_composite_external_potential.py b/python/tests/test_composite_external_potential.py index 9515921..e075944 100644 --- a/python/tests/test_composite_external_potential.py +++ b/python/tests/test_composite_external_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def comp(): diff --git a/python/tests/test_composite_flux.py b/python/tests/test_composite_flux.py index f1785eb..b92b7d0 100644 --- a/python/tests/test_composite_flux.py +++ b/python/tests/test_composite_flux.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def comp(): diff --git a/python/tests/test_composite_functional.py b/python/tests/test_composite_functional.py index ebd37f6..ec070a1 100644 --- a/python/tests/test_composite_functional.py +++ b/python/tests/test_composite_functional.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def comp(): diff --git a/python/tests/test_crank_nicolson_integrator.py b/python/tests/test_crank_nicolson_integrator.py index dfcba14..2a9153f 100644 --- a/python/tests/test_crank_nicolson_integrator.py +++ b/python/tests/test_crank_nicolson_integrator.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def cn(): diff --git a/python/tests/test_explicit_euler_integrator.py b/python/tests/test_explicit_euler_integrator.py index 2889542..1e0e3e6 100644 --- a/python/tests/test_explicit_euler_integrator.py +++ b/python/tests/test_explicit_euler_integrator.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def euler(): diff --git a/python/tests/test_exponential_wall_potential.py b/python/tests/test_exponential_wall_potential.py index 5da332b..18573b1 100644 --- a/python/tests/test_exponential_wall_potential.py +++ b/python/tests/test_exponential_wall_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def ew(): diff --git a/python/tests/test_external_field.py b/python/tests/test_external_field.py index b07c5d0..f896a48 100644 --- a/python/tests/test_external_field.py +++ b/python/tests/test_external_field.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def spherical_mesh_grand(): diff --git a/python/tests/test_field.py b/python/tests/test_field.py index 2409aff..c7f3d78 100644 --- a/python/tests/test_field.py +++ b/python/tests/test_field.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def field(): diff --git a/python/tests/test_grand_potential.py b/python/tests/test_grand_potential.py index 768cd2c..bbf9897 100644 --- a/python/tests/test_grand_potential.py +++ b/python/tests/test_grand_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - from .test_ideal_gas import f_ig, mu_ig from .test_rosenfeld_fmt import fex_py, muex_py diff --git a/python/tests/test_hard_wall_potential.py b/python/tests/test_hard_wall_potential.py index 697d82e..5730945 100644 --- a/python/tests/test_hard_wall_potential.py +++ b/python/tests/test_hard_wall_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def hw(): diff --git a/python/tests/test_harmonic_wall_potential.py b/python/tests/test_harmonic_wall_potential.py index 7d6dbd2..a700c2a 100644 --- a/python/tests/test_harmonic_wall_potential.py +++ b/python/tests/test_harmonic_wall_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def hw(): diff --git a/python/tests/test_implicit_euler_integrator.py b/python/tests/test_implicit_euler_integrator.py index aa39865..dc7bc7e 100644 --- a/python/tests/test_implicit_euler_integrator.py +++ b/python/tests/test_implicit_euler_integrator.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def euler(): diff --git a/python/tests/test_lennard_jones_93_wall_potential.py b/python/tests/test_lennard_jones_93_wall_potential.py index efa5692..ca70f83 100644 --- a/python/tests/test_lennard_jones_93_wall_potential.py +++ b/python/tests/test_lennard_jones_93_wall_potential.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def lj(): diff --git a/python/tests/test_mirror.py b/python/tests/test_mirror.py index f721114..1828282 100644 --- a/python/tests/test_mirror.py +++ b/python/tests/test_mirror.py @@ -1,6 +1,5 @@ -import pytest - import flyft.mirror +import pytest class _A: diff --git a/python/tests/test_parameter.py b/python/tests/test_parameter.py index 82631cc..3e91ea5 100644 --- a/python/tests/test_parameter.py +++ b/python/tests/test_parameter.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - class TimeParameter(flyft.CustomParameter): def __call__(self, state): diff --git a/python/tests/test_picard_iteration.py b/python/tests/test_picard_iteration.py index 5bc1cf5..b5ad3de 100644 --- a/python/tests/test_picard_iteration.py +++ b/python/tests/test_picard_iteration.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - from .test_ideal_gas import mu_ig from .test_rosenfeld_fmt import muex_py diff --git a/python/tests/test_rpy_diffusive_flux.py b/python/tests/test_rpy_diffusive_flux.py index 71d4524..89e433c 100644 --- a/python/tests/test_rpy_diffusive_flux.py +++ b/python/tests/test_rpy_diffusive_flux.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - @pytest.fixture def spherical_mesh_grand(): diff --git a/python/tests/test_state.py b/python/tests/test_state.py index 24924a6..a1aed37 100644 --- a/python/tests/test_state.py +++ b/python/tests/test_state.py @@ -1,6 +1,5 @@ -import pytest - import flyft +import pytest def test_init(state): diff --git a/python/tests/test_virial_expansion.py b/python/tests/test_virial_expansion.py index 4b5d0e1..3464b95 100644 --- a/python/tests/test_virial_expansion.py +++ b/python/tests/test_virial_expansion.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - def f_ex(B, rho): """Excess free-energy density of virial expansion""" diff --git a/python/tests/test_white_bear.py b/python/tests/test_white_bear.py index 2924ae2..6ef571b 100644 --- a/python/tests/test_white_bear.py +++ b/python/tests/test_white_bear.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - def fex_cs(eta, v): """Carnahan--Starling free-energy density of hard spheres""" diff --git a/python/tests/test_white_bear_mark_ii.py b/python/tests/test_white_bear_mark_ii.py index 73feed1..48ee0ec 100644 --- a/python/tests/test_white_bear_mark_ii.py +++ b/python/tests/test_white_bear_mark_ii.py @@ -1,8 +1,7 @@ +import flyft import numpy as np import pytest -import flyft - def fex_cs(eta, v): """Carnahan--Starling free-energy density of hard spheres""" diff --git a/python/validate/binary_mixture.py b/python/validate/binary_mixture.py index 4cdae6b..28e6620 100644 --- a/python/validate/binary_mixture.py +++ b/python/validate/binary_mixture.py @@ -1,6 +1,5 @@ -import numpy as np - import flyft +import numpy as np L = 50.0 dx = 0.2 diff --git a/python/validate/brownian_diffusion.py b/python/validate/brownian_diffusion.py index 8183677..0fc5d5a 100644 --- a/python/validate/brownian_diffusion.py +++ b/python/validate/brownian_diffusion.py @@ -1,6 +1,5 @@ -import numpy as np - import flyft +import numpy as np L = 10.0 dx = 0.05 diff --git a/python/validate/hard_sphere_evaporation.py b/python/validate/hard_sphere_evaporation.py index af01e9e..be85854 100644 --- a/python/validate/hard_sphere_evaporation.py +++ b/python/validate/hard_sphere_evaporation.py @@ -1,6 +1,5 @@ -import numpy as np - import flyft +import numpy as np L = 10.0 dx = 0.05 diff --git a/python/validate/hard_sphere_slit.py b/python/validate/hard_sphere_slit.py index 73d6fdd..bb2a014 100644 --- a/python/validate/hard_sphere_slit.py +++ b/python/validate/hard_sphere_slit.py @@ -1,6 +1,5 @@ -import numpy as np - import flyft +import numpy as np L = 16.0 diameters = {"A": 1.0} diff --git a/python/validate/stratification.py b/python/validate/stratification.py index a8cdb6f..ee04a00 100644 --- a/python/validate/stratification.py +++ b/python/validate/stratification.py @@ -1,6 +1,5 @@ -import numpy as np - import flyft +import numpy as np L = 38.5 dx = 0.05 diff --git a/src/data_layout.cc b/src/data_layout.cc index dc58eab..297809d 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -1,27 +1,37 @@ #include "flyft/data_layout.h" #include +#include namespace flyft { -DataLayout::DataLayout() : shape_(0) {} +DataLayout::DataLayout() {} -DataLayout::DataLayout(int shape) : shape_(shape) {} +DataLayout::DataLayout(const std::vector& shape) : shape_(shape) {} -int DataLayout::operator()(int idx) const +int DataLayout::operator()(const std::vector& idx) const { - return idx; + int N = idx.size(); + int one_d_idx = 0; + int stride = 1; + for (int k = N - 1; k >= 0; --k) + { + one_d_idx += idx[k] * stride; + stride *= shape_[k]; + } + + return one_d_idx; } -int DataLayout::shape() const +std::vector DataLayout::shape() const { return shape_; } int DataLayout::size() const { - return shape(); + return std::accumulate(shape_.begin(), shape_.end(), 1, std::multiplies()); } bool DataLayout::operator==(const DataLayout& other) const @@ -33,5 +43,4 @@ bool DataLayout::operator!=(const DataLayout& other) const { return !(*this == other); } - } // namespace flyft From 009f2d32a4705d5a03524b134940a0eaa5562ad7 Mon Sep 17 00:00:00 2001 From: Mayukh Kundu <79770572+mayukh33@users.noreply.github.com> Date: Thu, 12 Jun 2025 09:45:24 -0500 Subject: [PATCH 02/13] Mesh for anisotropic particles (#52) * Extend center, lower_bound, upper_bound, shape, step and bin to higher dimensions * Fix multi-dimensional indexing * Apply fix and add documentation * Minor fix * Apply review suggestions and minor fixes * Revert .pre-commit-config.yaml --- include/flyft/mesh.h | 153 +++++++++++++++++++++++++++++++++---------- src/mesh.cc | 119 +++++++++++++++++++++------------ 2 files changed, 198 insertions(+), 74 deletions(-) diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h index ab85cc0..6a24370 100644 --- a/include/flyft/mesh.h +++ b/include/flyft/mesh.h @@ -8,54 +8,111 @@ namespace flyft { - +//! Multidimensional mesh +/*! + * Multidimensional mesh maps N-dimensional indexing for + * the anistropic particles into one-dimensional flux. + */ class Mesh { public: + //! No default constructor Mesh() = delete; + + //! Constructor + /*! + * \param lower_bound Lower bound of the N-dimensional mesh. + * \param upper_bound Upper bound of the N-dimensional mesh. + * \param shape Shape of the index array. + * \param lower_bc Boundary condition of the lower bound. + * \param upper_bc Boundary condition of the upper bound. + */ Mesh(double lower_bound, double upper_bound, int shape, BoundaryType lower_bc, BoundaryType upper_bc); - std::shared_ptr slice(int start, int end) const; + //! Slice multidimensional index + /*! + * \param start Multidimensional start index. + * \param end Multidimensional end index. + * \return One-dimensional index. + * + */ + std::shared_ptr slice(const std::vector& start, const std::vector& end) const; //! Get position on the mesh, defined as center of bin - double center(int i) const; + /*! + * \param i Multidimensional bin index + * \return Center of the multidimensional bin + */ + std::vector center(const std::vector& i) const; //! Lower bound of entire mesh - double lower_bound() const; + /*! + * \return Lower bound of the multidimensional mesh + */ + std::vector lower_bound() const; //! Get lower bound of bin - double lower_bound(int i) const; + /*! + * \param i Multidimensional bin index + * \return Lower bound of the multidimensional bin + */ + std::vector lower_bound(const std::vector& i) const; //! Upper bound of entire mesh - double upper_bound() const; + /*! + * \return Upper bound of the multidimensional mesh + */ + std::vector upper_bound() const; //! Get upper bound of bin - double upper_bound(int i) const; + /*! + * \param i Multidimensional bin index + * \return Upper bound of the multidimensional bin + */ + std::vector upper_bound(const std::vector& i) const; //! 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 + /*! + * \return Volume of the mesh + */ virtual double volume() const = 0; //! Get volume of bin - virtual double volume(int i) const = 0; + /*! + * \param i Multidimensional bin index + * \return Volume of the bin + */ + virtual std::vector volume(int i) const = 0; //! Get the bin for a coordinate - int bin(double x) const; + /*! + * \param x Coordinate position of the bin + * \return Bin index + */ + std::vector bin(const std::vector& x) const; //! Length of the mesh - double L() const; + /*! + * \return Length of the mesh + */ + std::vector L() const; //! Shape of the mesh - int shape() const; + std::vector shape() const; //! Step size of the mesh - double step() const; + std::vector step() const; //! Boundary condition on lower bound of mesh BoundaryType lower_boundary_condition() const; @@ -63,35 +120,65 @@ class Mesh //! Boundary condition on upper bound of mesh 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; - - double integrateVolume(int idx, double f) const; - double integrateVolume(int idx, const DataView& f) const; - double integrateVolume(int idx, const DataView& f) const; - + //! Get the start index of the mesh + /*! + * \param dx Step size of the mesh + * \return Start index of the mesh + */ + std::vector asShape(const std::vector& dx) const; + + //! Get the length of the mesh + /*! + * \param shape Shape of the mesh + * \return Length of the mesh + */ + std::vector asLength(const std::vector& shape) const; + + //! 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; + + //! Get the integral of the volume over the mesh + /*! + * \param idx Multidimensional index + * \param f Function to integrate + * \return Integral of the function over the mesh + */ + double integrateVolume(const std::vector& idx, double f) const; + double integrateVolume(const std::vector& idx, const DataView& f) const; + double integrateVolume(const std::vector& idx, const DataView& f) const; + + //! 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; bool operator==(const Mesh& other) const; bool operator!=(const Mesh& other) const; protected: - double lower_; - int shape_; //!< Shape of the mesh - BoundaryType lower_bc_; - BoundaryType upper_bc_; - double step_; //!< Spacing between mesh points - int start_; + std::vector lower_; + std::vector shape_; //!< Shape of the mesh + BoundaryType lower_bc_; //!< Boundary condition of the lower bound + BoundaryType upper_bc_; //!< Boundary condition of the upper bound + std::vector step_; //!< Spacing between mesh points + std::vector start_; void validateBoundaryCondition() const; diff --git a/src/mesh.cc b/src/mesh.cc index 6e063f3..3303edc 100644 --- a/src/mesh.cc +++ b/src/mesh.cc @@ -4,19 +4,21 @@ namespace flyft { - -Mesh::Mesh(double lower_bound, - double upper_bound, - int shape, +Mesh::Mesh(std::vector lower_bound, + std::vector upper_bound, + std::vector shape, BoundaryType lower_bc, BoundaryType upper_bc) : lower_(lower_bound), shape_(shape), lower_bc_(lower_bc), upper_bc_(upper_bc), start_(0) { - step_ = (upper_bound - lower_bound) / shape_; + for (int idx = 0; idx < shape_.size(); ++idx) + { + step_.push_back((upper_bound[idx] - lower_bound[idx]) / shape_[idx]); + } validateBoundaryCondition(); } -std::shared_ptr Mesh::slice(int start, int end) const +std::shared_ptr Mesh::slice(const std::vector& start, const std::vector& end) const { if (lower_bc_ == BoundaryType::internal || upper_bc_ == BoundaryType::internal) { @@ -24,102 +26,137 @@ std::shared_ptr Mesh::slice(int start, int end) const } auto m = clone(); - m->start_ = start; - m->shape_ = end - start; - if (start > 0) - { - m->lower_bc_ = BoundaryType::internal; - } - if (end < shape_) + for (int idx = 0; idx < start.size(); ++idx) { - m->upper_bc_ = BoundaryType::internal; + m->start_[idx] = start[idx]; + m->shape_[idx] = end[idx] - start[idx]; + if (start[idx] > 0) + { + m->lower_bc_ = BoundaryType::internal; + } + if (end < shape_) + { + m->upper_bc_ = BoundaryType::internal; + } } + return m; } -double Mesh::center(int i) const +std::vector Mesh::center(const std::vector& i) const { - return lower_ + static_cast(start_ + i + 0.5) * step_; + std::vector temp; + for (int idx; idx < shape_.size(); ++idx) + { + temp.push_back(lower_[idx] + static_cast(start_[idx] + i[idx] + 0.5)); + } + return temp; } -int Mesh::bin(double x) const +std::vector Mesh::bin(const std::vector& x) const { - return static_cast((x - lower_) / step_) - start_; + std::vector temp; + for (int idx; idx < shape_.size(); ++idx) + { + temp.push_back(static_cast((x[idx] - lower_[idx]) / step_[idx]) - start_[idx]); + } + return temp; } -double Mesh::lower_bound() const +std::vector Mesh::lower_bound() const { - return lower_bound(0); + std::vector temp(shape_.size(), 0); + return lower_bound(temp); } -double Mesh::lower_bound(int i) const +std::vector Mesh::lower_bound(const std::vector& i) const { - return lower_ + (start_ + i) * step_; + std::vector temp; + for (int idx; idx < shape_.size(); ++idx) + { + temp.push_back(lower_[idx] + (start_[idx] + i[idx]) * step_[idx]); + } + return temp } -double Mesh::upper_bound() const +std::vector Mesh::upper_bound() const { return lower_bound(shape_); } -double Mesh::upper_bound(int i) const +std::vector Mesh::upper_bound(const std::vector& i) const { - return lower_bound(i + 1); + std::vector temp; + for (int idx; idx < shape_.size(); ++idx) + { + temp.push_back(i[idx] + 1); + } + return lower_bound(temp); } -double Mesh::L() const +std::vector Mesh::L() const { return asLength(shape_); } -int Mesh::shape() const +std::vector Mesh::shape() const { return shape_; } -double Mesh::step() const +std::vector Mesh::step() const { return step_; } -int Mesh::asShape(double dx) const +std::vector Mesh::asShape(const std::vector& dx) const { - return static_cast(std::ceil(dx / step_)); + std::vector temp; + for (int idx = 0; idx < shape_.size(); ++idx) + { + temp.push_back(static_cast(std::ceil(dx[idx] / step_[idx]))); + } + return temp; } -double Mesh::asLength(int shape) const +std::vector Mesh::asLength(const std::vector& shape) const { - return shape * step_; + std::vector temp; + for (int idx = 0; idx < shape_.size(); ++idx) + { + temp.push_back(shape[idx] * step_[idx]); + } + return temp; } -double Mesh::integrateSurface(int idx, double j_lo, double j_hi) const +double Mesh::integrateSurface(const std::vector& idx, double j_lo, double j_hi) const { return area(idx) * j_lo - area(idx + 1) * j_hi; } -double Mesh::integrateSurface(int idx, const DataView& j) const +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const { return integrateSurface(idx, j(idx), j(idx + 1)); } -double Mesh::integrateSurface(int idx, const DataView& j) const +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const { - return integrateSurface(idx, j(idx), j(idx + 1)); + return integrateSurface(const std::vector& idx, j(idx), j(idx + 1)); } -double Mesh::integrateVolume(int idx, double f) const +double Mesh::integrateVolume(std::vector idx, double f) const { return volume(idx) * f; } -double Mesh::integrateVolume(int idx, const DataView& f) const +double Mesh::integrateVolume(const std::vector& idx, const DataView& f) const { - return integrateVolume(idx, f(idx)); + return integrateVolume(const std::vector& idx, f(idx)); } -double Mesh::integrateVolume(int idx, const DataView& f) const +double Mesh::integrateVolume(const std::vector& idx, const DataView& f) const { - return integrateVolume(idx, f(idx)); + return integrateVolume(const std::vector& idx, f(idx)); } double Mesh::gradient(int idx, const DataView& f) const From 0784e0b0b7df1d968ecf0ba74e9af310427e1b83 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 12 Jun 2025 22:09:14 -0500 Subject: [PATCH 03/13] Fix build issues --- include/flyft/mesh.h | 2 ++ python/CMakeLists.txt | 1 + python/flyft/CMakeLists.txt | 6 ++---- python/tests/CMakeLists.txt | 6 ++---- src/mesh.cc | 2 ++ 5 files changed, 9 insertions(+), 8 deletions(-) diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h index 6a24370..0272dbc 100644 --- a/include/flyft/mesh.h +++ b/include/flyft/mesh.h @@ -33,6 +33,8 @@ class Mesh BoundaryType lower_bc, BoundaryType upper_bc); + virtual ~Mesh(); + //! Slice multidimensional index /*! * \param start Multidimensional start index. 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/mesh.cc b/src/mesh.cc index 3303edc..5782b04 100644 --- a/src/mesh.cc +++ b/src/mesh.cc @@ -18,6 +18,8 @@ Mesh::Mesh(std::vector lower_bound, validateBoundaryCondition(); } +Mesh::~Mesh() {} + std::shared_ptr Mesh::slice(const std::vector& start, const std::vector& end) const { if (lower_bc_ == BoundaryType::internal || upper_bc_ == BoundaryType::internal) From a77e8937f8eb179931f4797b90ec0c4b7f6d58f3 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 12 Jun 2025 22:10:07 -0500 Subject: [PATCH 04/13] Try to fix pre-commit warning --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1b63521..b91afcf 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -17,7 +17,7 @@ repos: - id: clang-format types_or: [c, c++] - repo: https://github.com/PyCQA/isort - rev: 5.13.2 + rev: 6.0.1 hooks: - id: isort - repo: https://github.com/psf/black-pre-commit-mirror From 15c4e4006499be7bc915cb2ab86591d26d0ed5f0 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 26 Jun 2025 11:37:59 -0500 Subject: [PATCH 05/13] Refactor data layout --- .vscode/extensions.json | 1 - include/flyft/data_layout.h | 92 ++++++++++++----- src/data_layout.cc | 196 +++++++++++++++++++++++++++++++++--- 3 files changed, 250 insertions(+), 39 deletions(-) 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/include/flyft/data_layout.h b/include/flyft/data_layout.h index d32077d..3787eee 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -1,55 +1,101 @@ #ifndef FLYFT_DATA_LAYOUT_H_ #define FLYFT_DATA_LAYOUT_H_ -#include -#include -#include -#include +#include namespace flyft { -//! Multidimensional array layout + +//! Memory layout of a multidimensional array. /*! - * DataLayout maps an N-dimensional index to a one-dimensional - * index using row-major ordering. Negative indexes are supported - * and interpreted as relative to the zero-index. + * 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 multi-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. */ class DataLayout { public: + //! Empty constructor DataLayout(); + //! Constructor /*! - * \param shape Shape of the index array. + * \param shape Shape of the array. + * \param num_dimensions Number of dimensions. */ - explicit DataLayout(const std::vector& shape); + DataLayout(const int* shape, char num_dimensions); + + //! Copy constructor + DataLayout(const DataLayout& other); + + //! Copy assignment operator + DataLayout& operator=(const DataLayout& other); + + //! Move constructor + DataLayout(const DataLayout&& other); - //! Map a multidimensional index to a one-dimensional index + //! Move assignment operator + DataLayout& operator=(const DataLayout&& other); + + //! Destructor + ~DataLayout(); + + //! Map a multidimensional index to a one-dimensional index. /*! - * \param idx Multidimensional index. + * \param multi_index Multidimensional index. * \return One-dimensional index. - * - * Row-major ordering is used. */ - int operator()(const std::vector& idx) const; + size_t operator()(const int* multi_index) const; - //! Shape of the layout in each dimension - std::vector shape() 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; - //! Total number of elements in the layout - int size() 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; - //! Test if two layouts are equal + //! Map a one-dimensional index to a multidimensional index offset. /*! - * \return True if the layouts have the same shape. + * \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; + + //! Shape of each dimension. + const int* shape() const; + + //! Number of dimensions. + char num_dimensions() const; + + //! Total number of elements. + size_t size() const; + + //! Test if two layouts are equal. bool operator==(const DataLayout& other) const; - //! Test if two layouts are not equal + //! Test if two layouts are not equal. bool operator!=(const DataLayout& other) const; private: - std::vector shape_; //!< Multidimensional shape of layout + int* shape_; //!< Shape of array along each dimension. + char num_dimensions_; //!< Number of dimensions. + size_t size_; //!< Total number of elements. + + //! Reset shape and dimensionality of layout + void reset(const int* shape, char num_dimensions); }; } // namespace flyft diff --git a/src/data_layout.cc b/src/data_layout.cc index 297809d..4987941 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -1,46 +1,212 @@ #include "flyft/data_layout.h" #include -#include namespace flyft { -DataLayout::DataLayout() {} +DataLayout::DataLayout() : shape_(nullptr), num_dimensions_(0), size_(0) {} -DataLayout::DataLayout(const std::vector& shape) : shape_(shape) {} +DataLayout::DataLayout(const int* shape, char num_dimensions) + { + reset(shape, num_dimensions); + } + +DataLayout::DataLayout(const DataLayout& other) + { + reset(other.shape_, other.num_dimensions_); + } + +DataLayout& DataLayout::operator=(const DataLayout& other) + { + if (this != &other) + { + reset(other.shape_, other.num_dimensions_); + } + return *this; + } + +DataLayout::DataLayout(const DataLayout&& other) + : shape_(std::move(other.shape_)), num_dimensions_(std::move(other.num_dimensions_)), + size_(std::move(other.size_)) + { + } + +DataLayout& DataLayout::operator=(const DataLayout&& other) + { + shape_ = std::move(other.shape_); + num_dimensions_ = std::move(other.num_dimensions_); + size_ = std::move(other.size_); + return *this; + } + +DataLayout::~DataLayout() + { + delete[] shape_; + } + +size_t DataLayout::operator()(const int* multi_index) const + { + size_t flat_index; + if (num_dimensions_ > 1) + { + // multiply out using row-major ordering + flat_index = 0; + size_t stride = 1; + for (char dim = num_dimensions_ - 1; dim >= 0; --dim) + { + flat_index += multi_index[dim] * stride; + stride *= shape_[dim]; + } + } + else + { + flat_index = multi_index[0]; + } + return flat_index; + } + +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 + size_t flat_index; + if (num_dimensions_ > 1) + { + // multiply out using row-major ordering + flat_index = 0; + size_t stride = 1; + for (char dim = num_dimensions_ - 1; dim >= 0; --dim) + { + flat_index += (offset[dim] + multi_index[dim]) * stride; + stride *= shape_[dim]; + } + } + else + { + flat_index = offset[0] + multi_index[0]; + } + return flat_index; + } -int DataLayout::operator()(const std::vector& idx) const +void DataLayout::operator()(int* multi_index, size_t flat_index) const { - int N = idx.size(); - int one_d_idx = 0; - int stride = 1; - for (int k = N - 1; k >= 0; --k) + if (num_dimensions_ > 1) { - one_d_idx += idx[k] * stride; - stride *= shape_[k]; + // repeatedly floor divide the flat index by stride associated with it + size_t idx = flat_index; + size_t stride = size_; + for (char dim = 0; dim < num_dimensions_ - 1; ++dim) + { + stride /= shape_[dim]; + multi_index[dim] = idx / stride; + idx -= multi_index[dim] * stride; + } } + else + { + multi_index[0] = flat_index; + } + } - return one_d_idx; +void DataLayout::operator()(int* multi_index, size_t flat_index, const int* offset) const + { + this->operator()(multi_index, flat_index); + if (offset) + { + for (char dim = 0; dim < num_dimensions_ - 1; ++dim) + { + multi_index[dim] -= offset[dim]; + } + } } -std::vector DataLayout::shape() const +const int* DataLayout::shape() const { return shape_; } -int DataLayout::size() const +char DataLayout::num_dimensions() const { - return std::accumulate(shape_.begin(), shape_.end(), 1, std::multiplies()); + return num_dimensions_; + } + +size_t DataLayout::size() const + { + return size_; } bool DataLayout::operator==(const DataLayout& other) const { - return (shape_ == other.shape_); + // dimensionality & size must match as first pass + bool same_layout = (num_dimensions_ == other.num_dimensions_ && size_ == other.size_); + if (same_layout) + { + // then, shape must match along all dimensions + if (shape_ && other.shape_) + { + for (char i = 0; i < num_dimensions_; ++i) + { + if (shape_[i] != other.shape_[i]) + { + same_layout = false; + } + } + } + // or, both pointers must be null + else if (!(!shape_ && !other.shape_)) + { + same_layout = false; + } + } + + return same_layout; } bool DataLayout::operator!=(const DataLayout& other) const { return !(*this == other); } + +void DataLayout::reset(const int* shape, char num_dimensions) + { + if (shape && num_dimensions > 0) + { + // free shape if it is allocated but the wrong size + if (shape_ && num_dimensions_ != num_dimensions) + { + delete shape_; + } + + // (re-)allocate shape memory if needed + if (!shape_) + { + shape_ = new int[num_dimensions]; + } + + // always copy shape values and recompute size + num_dimensions_ = num_dimensions; + size_ = 1; + for (char i = 0; i < num_dimensions; ++i) + { + shape_[i] = shape[i]; + size_ *= shape[i]; + } + } + else + { + if (shape_) + { + delete shape_; + } + num_dimensions_ = 0; + size_ = 0; + } + } + } // namespace flyft From ba84363f742b842b7b2f3a56c38c3e7414ef3eb6 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 26 Jun 2025 11:38:37 -0500 Subject: [PATCH 06/13] Start modifying data view --- include/flyft/data_view.h | 149 +++++++++++++++++--------------------- 1 file changed, 65 insertions(+), 84 deletions(-) diff --git a/include/flyft/data_view.h b/include/flyft/data_view.h index 4449d20..245e542 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -3,39 +3,43 @@ #include "flyft/data_layout.h" -#include -#include - namespace flyft { -template -//! Multidimensional array view +//! Accessor to a multidimensional array. /*! - * DataView provides ability to view a certain section of the multidimensional array + * A Dataview interprets a pointer as a multidimensional array according to a particular DataLayout. + * It provides random access to elements by multi-index, as well as a forward iterator. + * + * The view can be scoped to a subset of only some of the valid multi-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 + * multi-index of 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. */ +template class DataView { public: - using value_type = typename std::remove_reference::type; // Extracts datatype - using pointer = value_type*; // Pointer to the datatype - using reference = value_type&; // Reference to the datatype + using value_type = typename std::remove_reference::type; // start, std::vector end) + DataView(pointer data, const DataLayout& layout, const int* start, const int* end) : data_(data), layout_(layout), start_(start), end_(end) { } - //! Reference to the referenced element of the data. + // non-copyable, non-movable + DataView(const DataView& other) = delete; + DataView& operator=(const DataView& other) = delete; + DataView(const DataView&& other) = delete; + DataView& operator=(const DataView&& other) = delete; + + //! Destructor. + ~DataView() {} + + //! Access an element of the data at a given multidimensional index. /*! - * \param idx Array of indices for the multiple dimensions - * \return reference to the referenced data + * \param multi_index Multidimensional index to access. + * \return Data element. */ - reference operator()(const std::vector& idx) const + reference operator()(const int* multi_index) const { - std::vector temp; - std::transform(idx.begin(), idx.end(), start_.begin(), temp.begin(), std::plus()); - return data_[layout_(temp)]; + return data_[layout_(multi_index, start_)]; } + //! Number of dimensions. + char num_dimensions() const + { + return layout_.num_dimensions(); + } + +#if 0 //! Shape of the multi-dimensional array std::vector shape() const { - if (start_.size() != end_.size()) - { - throw std::invalid_argument( - "Arrays must have the same size for element-wise subtraction."); - } - - std::vector& result(start_.size()); - for (int i = 0; i < start_.size(); ++i) + std::vector& result(layout_.num_dimensions()); + for (char dim = 0; dim < layout_.num_dimensions(); ++dim) { - result[i] = start_[i] - end_[i]; + result[i] = end_[dim] - start_[dim]; } return result; } @@ -155,10 +151,11 @@ class DataView { return std::accumulate(shape().begin(), shape().end(), 1, std::multiplies()); } +#endif - //! Test if pointer to the data is not a null pointer + //! Test if view is not null /*! - * \return DataView is not null if it is a view of valid data. + * \returns True if data pointer being viewed is not null. */ explicit operator bool() const { @@ -183,27 +180,11 @@ class DataView return Iterator(*this, shape()); } - //! Test if two views are equal - /*! - * \return True if the view have the same pointer, layout, start and end points. - */ - bool operator==(const DataView& other) const - { - return (data_ == other.data_ && layout_ == other.layout_ && start_ == other.start_ - && end_ == other.end_); - } - - //! Test if two views are not equal - bool operator!=(const DataView& other) const - { - return !(*this == other); - } - private: - pointer data_; //!< Pointer to the array - DataLayout layout_; //!< Multidimensional array layout - std::vector start_; //!< Start of the array - std::vector end_; //!< End of the array + pointer data_; //!< Pointer to the array + DataLayout layout_; //!< Multidimensional array layout + int* start_; //!< Start of scoped view + int* end_; //!< End of scoped view }; } // namespace flyft From b71de6ea60a249a76fd7cb0800c3301656479693 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 26 Jun 2025 14:57:25 -0500 Subject: [PATCH 07/13] Draft updates to data view --- include/flyft/data_layout.h | 14 ++-- include/flyft/data_view.h | 142 +++++++++++++++++++++++++----------- src/data_layout.cc | 30 ++++---- 3 files changed, 124 insertions(+), 62 deletions(-) diff --git a/include/flyft/data_layout.h b/include/flyft/data_layout.h index 3787eee..83fb535 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -22,10 +22,10 @@ class DataLayout //! Constructor /*! - * \param shape Shape of the array. * \param num_dimensions Number of dimensions. + * \param shape Number of elements per dimension. */ - DataLayout(const int* shape, char num_dimensions); + DataLayout(char num_dimensions, const int* shape); //! Copy constructor DataLayout(const DataLayout& other); @@ -74,12 +74,12 @@ class DataLayout */ void operator()(int* multi_index, size_t flat_index, const int* offset) const; - //! Shape of each dimension. - const int* shape() const; - //! Number of dimensions. char num_dimensions() const; + //! Number of elements per dimension. + const int* shape() const; + //! Total number of elements. size_t size() const; @@ -90,12 +90,12 @@ class DataLayout bool operator!=(const DataLayout& other) const; private: - int* shape_; //!< Shape of array along each dimension. char num_dimensions_; //!< Number of dimensions. + int* shape_; //!< Number of elements per dimension. size_t size_; //!< Total number of elements. //! Reset shape and dimensionality of layout - void reset(const int* shape, char num_dimensions); + void reset(char num_dimensions, const int* shape); }; } // namespace flyft diff --git a/include/flyft/data_view.h b/include/flyft/data_view.h index 245e542..d4341df 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -8,7 +8,7 @@ namespace flyft //! Accessor to a multidimensional array. /*! - * A Dataview interprets a pointer as a multidimensional array according to a particular DataLayout. + * A DataView interprets a pointer as a multidimensional array according to a particular DataLayout. * It provides random access to elements by multi-index, as well as a forward iterator. * * The view can be scoped to a subset of only some of the valid multi-indexes. This is useful for @@ -36,10 +36,42 @@ class DataView //! Empty constructor Iterator() : Iterator(DataView()) {} - //! Construct from view - explicit Iterator(const DataView& view) : Iterator(view, 0) {} + //! Construct from view. + explicit Iterator(const DataView& view) : Iterator(view, nullptr) {} - Iterator(const DataView& view, size_t current) : view_(view), current_(current) {} + //! Construct from view starting from multidimensional index. + Iterator(const DataView& view, int* current) : view_(view) + { + current_ = new int[view_.num_dimensions()]; + for (char dim = 0; dim < view_.num_dimensions(); ++dim) + { + current_[dim] = (current) ? current[dim] : 0; + } + } + + // non-copyable + Iterator(const Iterator& other) = delete; + Iterator& operator=(const Iterator& other) = delete; + + //! Move constructor + Iterator(const Iterator&& other) + : view_(std::move(other.view_)), current_(std::move(other.current_)) + { + } + + //! Move assignment operator. + Iterator& operator=(const Iterator&& other) + { + view_ = std::move(other.view_); + current_ = std::move(other.current_); + return *this; + } + + //! Destructor. + ~Iterator() + { + delete[] current_; + } reference operator*() const { @@ -58,15 +90,19 @@ class DataView Iterator& operator++() { - ++current_; - return *this; - } + // increment the multi-index of the last dimension by 1 + char dim = view_.num_dimensions() - 1; + ++current_[dim]; + + // then, carry forward + const auto view_shape = view_.shape(); + while (current_[dim] == view_shape[dim] && dim >= 1) + { + current_[dim] = 0; + ++current_[dim--]; + } - Iterator operator++(int) - { - Iterator tmp(*this); - ++current_; - return tmp; + return *this; } bool operator==(const Iterator& other) const @@ -81,20 +117,27 @@ class DataView private: DataView view_; - size_t current_flat_index; + int* current_; }; //! Empty constructor. - DataView() : data_(nullptr), layout_(DataLayout()), start_(nullptr), end_(nullptr) {} + DataView() : data_(nullptr), layout_(DataLayout()), start_(nullptr), shape_(nullptr) {} //! Constructor. /*! * \param data Pointer to the data. * \param layout Layout of the data. */ - DataView(pointer data, const DataLayout& layout) - : data_(data), layout_(layout), start_(nullptr), end_(nullptr) + DataView(pointer data, const DataLayout& layout) : data_(data), layout_(layout) { + start_ = new int[layout_.num_dimensions()]; + shape_ = new int[layout_.num_dimensions()]; + const auto layout_shape = layout._shape(); + for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + { + start_[dim] = 0; + shape_[dim] = layout_shape[dim]; + } } //! Scoped-view constructor. @@ -102,11 +145,18 @@ class DataView * \param data_ Pointer to the data. * \param layout_ Layout of the data. * \param start_ Lower bound of multidimensional indexes within scope. - * \param end_ Upper bound of multidimensional indexes within scope. + * \param shape_ Shape of multidimensional indexes within scope, beginning with start. */ - DataView(pointer data, const DataLayout& layout, const int* start, const int* end) - : data_(data), layout_(layout), start_(start), end_(end) + DataView(pointer data, const DataLayout& layout, const int* start, const int* shape) + : data_(data), layout_(layout) { + start_ = new int[layout_.num_dimensions()]; + shape_ = new int[layout_.num_dimensions()]; + for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + { + start_[dim] = start[dim]; + shape_[dim] = shape[dim]; + } } // non-copyable, non-movable @@ -116,7 +166,11 @@ class DataView DataView& operator=(const DataView&& other) = delete; //! Destructor. - ~DataView() {} + ~DataView() + { + delete[] start_; + delete[] shape_; + } //! Access an element of the data at a given multidimensional index. /*! @@ -134,24 +188,30 @@ class DataView return layout_.num_dimensions(); } -#if 0 - //! Shape of the multi-dimensional array - std::vector shape() const + //! Number of elements per dimension in the scope of the view. + const int* shape() const { - std::vector& result(layout_.num_dimensions()); - for (char dim = 0; dim < layout_.num_dimensions(); ++dim) - { - result[i] = end_[dim] - start_[dim]; - } - return result; + return shape_; } - //! Number of elements in the multi-dimensional array - int size() const + //! Number of elements in the scope of the view. + size_t size() const { - return std::accumulate(shape().begin(), shape().end(), 1, std::multiplies()); + size_t size; + if (shape_) + { + size = 1; + for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + { + size *= shape_[dim]; + } + } + else + { + size = 0; + } + return size; } -#endif //! Test if view is not null /*! @@ -164,27 +224,27 @@ class DataView //! Start point of the data /*! - * \return data at index position 0. + * \return Iterator to start of data. */ Iterator begin() const { - return Iterator(*this, 0); + return Iterator(*this); } //! End point of the data /*! - * \return data at end of the array. + * \return Iterator to end of data. */ Iterator end() const { - return Iterator(*this, shape()); + return Iterator(*this, shape_); } private: - pointer data_; //!< Pointer to the array - DataLayout layout_; //!< Multidimensional array layout - int* start_; //!< Start of scoped view - int* end_; //!< End of scoped view + pointer data_; //!< Data. + DataLayout layout_; //!< Multidimensional array layout. + int* start_; //!< Multi-index for start of scoped view. + int* shape_; //!< Number of elements per dimension of scoped view. }; } // namespace flyft diff --git a/src/data_layout.cc b/src/data_layout.cc index 4987941..59c7c6b 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -5,37 +5,37 @@ namespace flyft { -DataLayout::DataLayout() : shape_(nullptr), num_dimensions_(0), size_(0) {} +DataLayout::DataLayout() : num_dimensions_(0), shape_(nullptr), size_(0) {} -DataLayout::DataLayout(const int* shape, char num_dimensions) +DataLayout::DataLayout(char num_dimensions, const int* shape) { - reset(shape, num_dimensions); + reset(num_dimensions, shape); } DataLayout::DataLayout(const DataLayout& other) { - reset(other.shape_, other.num_dimensions_); + reset(other.num_dimensions_, other.shape_); } DataLayout& DataLayout::operator=(const DataLayout& other) { if (this != &other) { - reset(other.shape_, other.num_dimensions_); + reset(other.num_dimensions_, other.shape_); } return *this; } DataLayout::DataLayout(const DataLayout&& other) - : shape_(std::move(other.shape_)), num_dimensions_(std::move(other.num_dimensions_)), + : num_dimensions_(std::move(other.num_dimensions_)), shape_(std::move(other.shape_)), size_(std::move(other.size_)) { } DataLayout& DataLayout::operator=(const DataLayout&& other) { - shape_ = std::move(other.shape_); num_dimensions_ = std::move(other.num_dimensions_); + shape_ = std::move(other.shape_); size_ = std::move(other.size_); return *this; } @@ -126,14 +126,14 @@ void DataLayout::operator()(int* multi_index, size_t flat_index, const int* offs } } -const int* DataLayout::shape() const +char DataLayout::num_dimensions() const { - return shape_; + return num_dimensions_; } -char DataLayout::num_dimensions() const +const int* DataLayout::shape() const { - return num_dimensions_; + return shape_; } size_t DataLayout::size() const @@ -173,14 +173,15 @@ bool DataLayout::operator!=(const DataLayout& other) const return !(*this == other); } -void DataLayout::reset(const int* shape, char num_dimensions) +void DataLayout::reset(char num_dimensions, const int* shape) { if (shape && num_dimensions > 0) { // free shape if it is allocated but the wrong size if (shape_ && num_dimensions_ != num_dimensions) { - delete shape_; + delete[] shape_; + shape_ = nullptr; } // (re-)allocate shape memory if needed @@ -202,7 +203,8 @@ void DataLayout::reset(const int* shape, char num_dimensions) { if (shape_) { - delete shape_; + delete[] shape_; + shape_ = nullptr; } num_dimensions_ = 0; size_ = 0; From cf8df039e7e921d669c3332e9194adc275c4d8bb Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 26 Jun 2025 14:58:04 -0500 Subject: [PATCH 08/13] Draft updates to field --- include/flyft/field.h | 108 ++++++++++++++++++++++++++++-------------- 1 file changed, 73 insertions(+), 35 deletions(-) diff --git a/include/flyft/field.h b/include/flyft/field.h index a32e974..9e77dbb 100644 --- a/include/flyft/field.h +++ b/include/flyft/field.h @@ -21,11 +21,15 @@ class GenericField : public TrackedObject 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) + + explicit GenericField(char num_dimensions, const int* shape) + : GenericField(num_dimensions, shape, nullptr) { - reshape(shape, buffer_shape); + } + + GenericField(char num_dimensions, const int* shape, const int* buffer_shape) + { + reshape(num_dimensions, shape, buffer_shape); } // noncopyable / nonmovable @@ -36,70 +40,99 @@ class GenericField : public TrackedObject ~GenericField() { - if (data_) - delete[] data_; + delete[] data_; + delete[] shape_; + delete[] buffer_shape_; + delete[] full_shape_; } - T& operator()(int idx) + T& operator()(const int* multi_index) { - return view()(idx); + return view()(multi_index); } - const T& operator()(int idx) const + const T& operator()(const int* multi_index) const { - return const_view()(idx); + return const_view()(multi_index); } - int shape() const + const int* shape() const { return shape_; } - int buffer_shape() const + const int* buffer_shape() const { return buffer_shape_; } - int full_shape() const + const int* full_shape() const { - return shape_ + 2 * buffer_shape_; + return layout_.shape(); } View view() { token_.stageAndCommit(); - return View(data_, layout_, buffer_shape_, full_shape() - buffer_shape_); + return View(data_, layout_, buffer_shape_, shape_); } ConstantView const_view() const { - return ConstantView(static_cast(data_), - layout_, - buffer_shape_, - full_shape() - buffer_shape_); + return ConstantView(static_cast(data_), layout_, buffer_shape_, shape_); } View full_view() { token_.stageAndCommit(); - return View(data_, layout_, 0, full_shape()); + return View(data_, layout_); } ConstantView const_full_view() const { - return ConstantView(data_, layout_, 0, full_shape()); + return ConstantView(data_, layout_); } - void reshape(int shape, int buffer_shape) + void reshape(char num_dimensions, const int* shape, const int* buffer_shape) { - if (shape != shape_ || buffer_shape != buffer_shape_) + if (num_dimensions != layout_.num_dimensions()) + { + delete[] shape_; + delete[] buffer_shape_; + delete[] full_shape_; + + shape_ = new int[num_dimensions] {}; + buffer_shape_ = new int[num_dimensions] {}; + full_shape_ = new int[num_dimensions] {}; + } + + bool same_shape = true; + bool same_buffer_shape = true; + for (char dim = 0; dim < num_dimensions; ++dim) + { + if (shape[dim] != shape_[dim]) + { + same_shape = false; + shape_[dim] = shape[dim]; + } + + if (buffer_shape[dim] != buffer_shape_[dim]) + { + same_buffer_shape = false; + buffer_shape_[dim] = buffer_shape[dim]; + } + + full_shape_[dim] = shape_[dim] + 2 * buffer_shape_[dim]; + } + + if (!(same_shape && same_buffer_shape)) { - DataLayout layout(shape + 2 * buffer_shape); + DataLayout layout(num_dimensions, full_shape_); // if data is alloc'd, decide what to do with it - if (data_ != nullptr) + if (data_) { - if (shape == shape_) + if (same_shape) { // data shape is the same but buffer changed, copy T* tmp = new T[layout.size()]; @@ -107,7 +140,7 @@ class GenericField : public TrackedObject // copy the data from old to new auto view_old = const_view(); - DataView view_new(tmp, layout, buffer_shape, buffer_shape + shape); + DataView view_new(tmp, layout, buffer_shape, shape); std::copy(view_old.begin(), view_old.end(), view_new.begin()); std::swap(data_, tmp); delete[] tmp; @@ -126,26 +159,30 @@ class GenericField : public TrackedObject } // if data is still not alloc'd, make it so - if (data_ == nullptr && layout.size() > 0) + if (!data_ && 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) + void setBuffer(const int* buffer_shape) { reshape(shape_, buffer_shape); } - void requestBuffer(int buffer_shape) + void requestBuffer(const int* buffer_shape) { - if (buffer_shape > buffer_shape_) + bool any_increase = false; + for (char dim = 0; !any_increase && dim < layout_.num_dimensions(); ++dim) + { + any_increase |= (buffer_shape[dim] > buffer_shape_[dim]); + } + + if (any_increase) { setBuffer(buffer_shape); } @@ -153,8 +190,9 @@ class GenericField : public TrackedObject private: T* data_; - int shape_; - int buffer_shape_; + int* shape_; + int* buffer_shape_; + int* full_shape_; DataLayout layout_; }; From fe743363e48c591aa2678d2b80692dd24b479bac Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Fri, 4 Jul 2025 12:12:58 -0500 Subject: [PATCH 09/13] Fix compile issues --- include/flyft/data_layout.h | 12 ++++++------ include/flyft/data_view.h | 14 ++++++++------ src/data_layout.cc | 18 +++++++++--------- 3 files changed, 23 insertions(+), 21 deletions(-) diff --git a/include/flyft/data_layout.h b/include/flyft/data_layout.h index 83fb535..f5c3ce2 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -25,7 +25,7 @@ class DataLayout * \param num_dimensions Number of dimensions. * \param shape Number of elements per dimension. */ - DataLayout(char num_dimensions, const int* shape); + DataLayout(int num_dimensions, const int* shape); //! Copy constructor DataLayout(const DataLayout& other); @@ -75,7 +75,7 @@ class DataLayout void operator()(int* multi_index, size_t flat_index, const int* offset) const; //! Number of dimensions. - char num_dimensions() const; + int num_dimensions() const; //! Number of elements per dimension. const int* shape() const; @@ -90,12 +90,12 @@ class DataLayout bool operator!=(const DataLayout& other) const; private: - char num_dimensions_; //!< Number of dimensions. - int* shape_; //!< Number of elements per dimension. - size_t size_; //!< Total number of elements. + int num_dimensions_; //!< Number of dimensions. + int* shape_; //!< Number of elements per dimension. + size_t size_; //!< Total number of elements. //! Reset shape and dimensionality of layout - void reset(char num_dimensions, const int* shape); + void reset(int num_dimensions, const int* shape); }; } // namespace flyft diff --git a/include/flyft/data_view.h b/include/flyft/data_view.h index d4341df..7f470fe 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -3,6 +3,8 @@ #include "flyft/data_layout.h" +#include + namespace flyft { @@ -43,7 +45,7 @@ class DataView Iterator(const DataView& view, int* current) : view_(view) { current_ = new int[view_.num_dimensions()]; - for (char dim = 0; dim < view_.num_dimensions(); ++dim) + for (int dim = 0; dim < view_.num_dimensions(); ++dim) { current_[dim] = (current) ? current[dim] : 0; } @@ -91,7 +93,7 @@ class DataView Iterator& operator++() { // increment the multi-index of the last dimension by 1 - char dim = view_.num_dimensions() - 1; + int dim = view_.num_dimensions() - 1; ++current_[dim]; // then, carry forward @@ -133,7 +135,7 @@ class DataView start_ = new int[layout_.num_dimensions()]; shape_ = new int[layout_.num_dimensions()]; const auto layout_shape = layout._shape(); - for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { start_[dim] = 0; shape_[dim] = layout_shape[dim]; @@ -152,7 +154,7 @@ class DataView { start_ = new int[layout_.num_dimensions()]; shape_ = new int[layout_.num_dimensions()]; - for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { start_[dim] = start[dim]; shape_[dim] = shape[dim]; @@ -183,7 +185,7 @@ class DataView } //! Number of dimensions. - char num_dimensions() const + int num_dimensions() const { return layout_.num_dimensions(); } @@ -201,7 +203,7 @@ class DataView if (shape_) { size = 1; - for (char dim = 0; dim < layout_.num_dimensions(); ++dim) + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { size *= shape_[dim]; } diff --git a/src/data_layout.cc b/src/data_layout.cc index 59c7c6b..56db9a9 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -7,7 +7,7 @@ namespace flyft DataLayout::DataLayout() : num_dimensions_(0), shape_(nullptr), size_(0) {} -DataLayout::DataLayout(char num_dimensions, const int* shape) +DataLayout::DataLayout(int num_dimensions, const int* shape) { reset(num_dimensions, shape); } @@ -53,7 +53,7 @@ size_t DataLayout::operator()(const int* multi_index) const // multiply out using row-major ordering flat_index = 0; size_t stride = 1; - for (char dim = num_dimensions_ - 1; dim >= 0; --dim) + for (int dim = num_dimensions_ - 1; dim >= 0; --dim) { flat_index += multi_index[dim] * stride; stride *= shape_[dim]; @@ -81,7 +81,7 @@ size_t DataLayout::operator()(const int* multi_index, const int* offset) const // multiply out using row-major ordering flat_index = 0; size_t stride = 1; - for (char dim = num_dimensions_ - 1; dim >= 0; --dim) + for (int dim = num_dimensions_ - 1; dim >= 0; --dim) { flat_index += (offset[dim] + multi_index[dim]) * stride; stride *= shape_[dim]; @@ -101,7 +101,7 @@ void DataLayout::operator()(int* multi_index, size_t flat_index) const // repeatedly floor divide the flat index by stride associated with it size_t idx = flat_index; size_t stride = size_; - for (char dim = 0; dim < num_dimensions_ - 1; ++dim) + for (int dim = 0; dim < num_dimensions_ - 1; ++dim) { stride /= shape_[dim]; multi_index[dim] = idx / stride; @@ -119,14 +119,14 @@ void DataLayout::operator()(int* multi_index, size_t flat_index, const int* offs this->operator()(multi_index, flat_index); if (offset) { - for (char dim = 0; dim < num_dimensions_ - 1; ++dim) + for (int dim = 0; dim < num_dimensions_ - 1; ++dim) { multi_index[dim] -= offset[dim]; } } } -char DataLayout::num_dimensions() const +int DataLayout::num_dimensions() const { return num_dimensions_; } @@ -150,7 +150,7 @@ bool DataLayout::operator==(const DataLayout& other) const // then, shape must match along all dimensions if (shape_ && other.shape_) { - for (char i = 0; i < num_dimensions_; ++i) + for (int i = 0; i < num_dimensions_; ++i) { if (shape_[i] != other.shape_[i]) { @@ -173,7 +173,7 @@ bool DataLayout::operator!=(const DataLayout& other) const return !(*this == other); } -void DataLayout::reset(char num_dimensions, const int* shape) +void DataLayout::reset(int num_dimensions, const int* shape) { if (shape && num_dimensions > 0) { @@ -193,7 +193,7 @@ void DataLayout::reset(char num_dimensions, const int* shape) // always copy shape values and recompute size num_dimensions_ = num_dimensions; size_ = 1; - for (char i = 0; i < num_dimensions; ++i) + for (int i = 0; i < num_dimensions; ++i) { shape_[i] = shape[i]; size_ *= shape[i]; From 78af7dc6cddbfe48451d4a912d72cdd3eb5f55cc Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Sat, 5 Jul 2025 12:36:26 -0500 Subject: [PATCH 10/13] First draft of changes to Mesh --- .vscode/settings.json | 2 +- include/flyft/cartesian_mesh.h | 35 +++- include/flyft/data_view.h | 3 +- include/flyft/mesh.h | 301 +++++++++++++++++++-------- include/flyft/spherical_mesh.h | 14 +- src/cartesian_mesh.cc | 50 ++++- src/mesh.cc | 362 ++++++++++++++++++++++++--------- src/spherical_mesh.cc | 44 ++-- 8 files changed, 578 insertions(+), 233 deletions(-) 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/data_view.h b/include/flyft/data_view.h index 7f470fe..4926c25 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -4,6 +4,7 @@ #include "flyft/data_layout.h" #include +#include namespace flyft { @@ -134,7 +135,7 @@ class DataView { start_ = new int[layout_.num_dimensions()]; shape_ = new int[layout_.num_dimensions()]; - const auto layout_shape = layout._shape(); + const auto layout_shape = layout_.shape(); for (int dim = 0; dim < layout_.num_dimensions(); ++dim) { start_[dim] = 0; diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h index 0272dbc..a2f1a9b 100644 --- a/include/flyft/mesh.h +++ b/include/flyft/mesh.h @@ -8,134 +8,231 @@ namespace flyft { -//! Multidimensional mesh +//! Mesh /*! - * Multidimensional mesh maps N-dimensional indexing for - * the anistropic particles into one-dimensional flux. + * 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 + // No default constructor Mesh() = delete; //! Constructor /*! - * \param lower_bound Lower bound of the N-dimensional mesh. - * \param upper_bound Upper bound of the N-dimensional mesh. - * \param shape Shape of the index array. - * \param lower_bc Boundary condition of the lower bound. - * \param upper_bc Boundary condition of the upper bound. + * \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(double lower_bound, + Mesh(int num_orientation_dim, + int* shape, + double lower_bound, double upper_bound, - int shape, BoundaryType lower_bc, BoundaryType upper_bc); + //! Copy constructor. + Mesh(const Mesh& other); + + //! Copy assignment operator. + Mesh& operator=(const Mesh& other); + + //! Move constructor. + Mesh(const Mesh&& other); + + //! Move assignment operator. + Mesh& operator=(const Mesh&& other); + virtual ~Mesh(); - //! Slice multidimensional index + //! Slice mesh by position coordinate. /*! - * \param start Multidimensional start index. - * \param end Multidimensional end index. - * \return One-dimensional index. + * \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(const std::vector& start, const std::vector& end) const; + std::shared_ptr slice(int start, int end) const; - //! Get position on the mesh, defined as center of bin + //! Calculate lower bound of cell. /*! - * \param i Multidimensional bin index - * \return Center of the multidimensional bin + * \param coordinate Lower bound of cell. + * \param cell Cell multi-index. */ - std::vector center(const std::vector& i) const; + void cell_lower_bound(double* coordinate, const int* cell) const; - //! Lower bound of entire mesh + //! Calculate lower bound of cell in a given dimension. /*! - * \return Lower bound of the multidimensional mesh + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns Lower bound of cell in specified dimension. */ - std::vector lower_bound() const; + double cell_lower_bound(int dimension, int index) const; - //! Get lower bound of bin + //! Calculate lower bound of cell in position space. /*! - * \param i Multidimensional bin index - * \return Lower bound of the multidimensional bin + * \param coordinate Lower bound of cell in position space. + * \param cell Cell multi-index. */ - std::vector lower_bound(const std::vector& i) const; + void cell_position_lower_bound(double* coordinate, const int* cell) const; - //! Upper bound of entire mesh + //! Calculate lower bound of cell in orientation space. /*! - * \return Upper bound of the multidimensional mesh + * \param coordinate Lower bound of cell in orientation space. + * \param cell Cell multi-index. */ - std::vector upper_bound() const; + void cell_orientation_lower_bound(double* coordinate, const int* cell) const; - //! Get upper bound of bin + //! Calculate center of cell. /*! - * \param i Multidimensional bin index - * \return Upper bound of the multidimensional bin + * \param coordinate Center of cell. + * \param cell Cell multi-index. */ - std::vector upper_bound(const std::vector& i) const; + void cell_center(double* coordinate, const int* cell) const; - //! Get surface area of lower edge of bin + //! Calculate center of cell in a given dimension. /*! - * \param i Multidimensional bin index - * \return Upper bound of the multidimensional bin + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns Center of cell in specified dimension. */ - virtual double area(int i) const = 0; + double cell_center(int dimension, int index) const; - //! Get volume of mesh + //! Calculate center of cell in position space. /*! - * \return Volume of the mesh + * \param coordinate Center of cell in position space. + * \param cell Cell multi-index. */ - virtual double volume() const = 0; + void cell_position_center(double* coordinate, const int* cell) const; - //! Get volume of bin + //! Calculate center of cell in orientation space. /*! - * \param i Multidimensional bin index - * \return Volume of the bin + * \param coordinate Center of cell in orientation space. + * \param cell Cell multi-index. */ - virtual std::vector volume(int i) const = 0; + void cell_orientation_center(double* coordinate, const int* cell) const; - //! Get the bin for a coordinate + //! Calculate upper bound of cell. /*! - * \param x Coordinate position of the bin - * \return Bin index + * \param coordinate Upper bound of cell. + * \param cell Cell multi-index. */ - std::vector bin(const std::vector& x) const; + void cell_upper_bound(double* coordinate, const int* cell) const; - //! Length of the mesh + //! Calculate upper bound of cell in a given dimension. /*! - * \return Length of the mesh + * \param dimension Dimension of coordinate. + * \param index Cell index in given dimension. + * \returns upper bound of cell in specified dimension. */ - std::vector L() const; + double cell_upper_bound(int dimension, int index) const; - //! Shape of the mesh - std::vector shape() 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; - //! Step size of the mesh - std::vector step() 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; - //! Boundary condition on lower bound of mesh - BoundaryType lower_boundary_condition() 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; +#endif - //! Boundary condition on upper bound of mesh - BoundaryType upper_boundary_condition() const; + //! 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; - //! Get the start index of the mesh + //! Calculate volume of cell in position space. /*! - * \param dx Step size of the mesh - * \return Start index of the mesh + * \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. */ - std::vector asShape(const std::vector& dx) const; + virtual double cell_position_volume(const int* cell) const = 0; - //! Get the length of the mesh + //! Calculate volume of cell in orientation space. /*! - * \param shape Shape of the mesh - * \return Length of the mesh + * \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. */ - std::vector asLength(const std::vector& shape) const; + 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 at position coordinate upper bound. + BoundaryType upper_boundary_condition() const; + +#if 0 //! Get the integral of the surface over the mesh /*! * \param shape Shape of the mesh @@ -144,17 +241,37 @@ class 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); + } - //! Get the integral of the volume over the mesh + //! Integrate function over cell orientation volume. /*! - * \param idx Multidimensional index - * \param f Function to integrate - * \return Integral of the function over the mesh + * \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. */ - double integrateVolume(const std::vector& idx, double f) const; - double integrateVolume(const std::vector& idx, const DataView& f) const; - double integrateVolume(const std::vector& idx, const DataView& f) const; + 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 @@ -170,23 +287,42 @@ class 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: - std::vector lower_; - std::vector shape_; //!< Shape of the mesh - BoundaryType lower_bc_; //!< Boundary condition of the lower bound - BoundaryType upper_bc_; //!< Boundary condition of the upper bound - std::vector step_; //!< Spacing between mesh points - std::vector 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 { @@ -211,6 +347,7 @@ typename std::remove_const::type Mesh::interpolate(double x, const DataView clone() const override; }; } // namespace flyft 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/mesh.cc b/src/mesh.cc index 5782b04..053f632 100644 --- a/src/mesh.cc +++ b/src/mesh.cc @@ -1,164 +1,334 @@ #include "flyft/mesh.h" +#include #include namespace flyft { -Mesh::Mesh(std::vector lower_bound, - std::vector upper_bound, - std::vector shape, +Mesh::Mesh(int num_orientation_dim, + int* shape, + double lower_bound, + double upper_bound, 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) { - for (int idx = 0; idx < shape_.size(); ++idx) + 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) { - step_.push_back((upper_bound[idx] - lower_bound[idx]) / shape_[idx]); + 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]; } - validateBoundaryCondition(); } -Mesh::~Mesh() {} - -std::shared_ptr Mesh::slice(const std::vector& start, const std::vector& 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]; } + } - auto m = clone(); - for (int idx = 0; idx < start.size(); ++idx) +Mesh& Mesh::operator=(const Mesh& other) + { + if (this != &other) { - m->start_[idx] = start[idx]; - m->shape_[idx] = end[idx] - start[idx]; - if (start[idx] > 0) - { - m->lower_bc_ = BoundaryType::internal; - } - if (end < shape_) + 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) { - m->upper_bc_ = BoundaryType::internal; + 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_[0] = end - start; + if (start > 0) + { + m->lower_bc_ = BoundaryType::internal; + } + if (end < shape_[0]) + { + m->upper_bc_ = BoundaryType::internal; + } return m; } -std::vector Mesh::center(const std::vector& i) const +void Mesh::cell_lower_bound(double* coordinate, const int* cell) const + { + cell_position_lower_bound(coordinate, cell); + cell_orientation_lower_bound(coordinate + num_position_dim_, cell); + } + +double Mesh::cell_lower_bound(int dimension, int index) const { - std::vector temp; - for (int idx; idx < shape_.size(); ++idx) + if (dimension < num_position_dim_) + { + assert(num_position_dim_ == 1); + return lower_[0] + (start_ + index) * step_[0]; + } + else { - temp.push_back(lower_[idx] + static_cast(start_[idx] + i[idx] + 0.5)); + return lower_[dimension] + index * step_[dimension]; } - return temp; } -std::vector Mesh::bin(const std::vector& x) const +void Mesh::cell_position_lower_bound(double* coordinate, const int* cell) const + { + assert(num_position_dim_ == 1); + coordinate[0] = cell_lower_bound(0, cell[0]); + } + +void Mesh::cell_orientation_lower_bound(double* coordinate, const int* cell) const { - std::vector temp; - for (int idx; idx < shape_.size(); ++idx) + for (int orientation_dim = 0; orientation_dim < num_orientation_dim_; ++orientation_dim) { - temp.push_back(static_cast((x[idx] - lower_[idx]) / step_[idx]) - start_[idx]); + coordinate[orientation_dim] + = cell_lower_bound(num_position_dim_ + orientation_dim, cell[orientation_dim]); } - return temp; } -std::vector Mesh::lower_bound() const +void Mesh::cell_center(double* coordinate, const int* cell) const { - std::vector temp(shape_.size(), 0); - return lower_bound(temp); + cell_position_center(coordinate, cell); + cell_orientation_center(coordinate + num_position_dim_, cell); } -std::vector Mesh::lower_bound(const std::vector& i) const +double Mesh::cell_center(int dimension, int index) const { - std::vector temp; - for (int idx; idx < shape_.size(); ++idx) + if (dimension < num_position_dim_) { - temp.push_back(lower_[idx] + (start_[idx] + i[idx]) * step_[idx]); + assert(num_position_dim_ == 1); + return lower_[0] + (start_ + index + 0.5) * step_[0]; + } + else + { + return lower_[dimension] + (index + 0.5) * step_[dimension]; } - return temp } -std::vector Mesh::upper_bound() const +void Mesh::cell_position_center(double* coordinate, const int* cell) const { - return lower_bound(shape_); + assert(num_position_dim_ == 1); + coordinate[0] = cell_center(0, cell[0]); } -std::vector Mesh::upper_bound(const std::vector& i) const +void Mesh::cell_orientation_center(double* coordinate, const int* cell) const { - std::vector temp; - for (int idx; idx < shape_.size(); ++idx) + for (int orientation_dim = 0; orientation_dim < num_orientation_dim_; ++orientation_dim) { - temp.push_back(i[idx] + 1); + coordinate[orientation_dim] + = cell_center(num_position_dim_ + orientation_dim, cell[orientation_dim]); } - return lower_bound(temp); } -std::vector Mesh::L() const +void Mesh::cell_upper_bound(double* coordinate, const int* cell) const { - return asLength(shape_); + cell_position_upper_bound(coordinate, cell); + cell_orientation_upper_bound(coordinate + num_position_dim_, cell); } -std::vector Mesh::shape() const +double Mesh::cell_upper_bound(int dimension, int index) const { - return shape_; + 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]; + } } -std::vector Mesh::step() const +void Mesh::cell_position_upper_bound(double* coordinate, const int* cell) const { - return step_; + assert(num_position_dim_ == 1); + coordinate[0] = cell_upper_bound(0, cell[0]); + } + +void Mesh::cell_orientation_upper_bound(double* coordinate, const int* cell) const + { + 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::cell_volume(const int* cell) const + { + return cell_position_volume(cell) * cell_orientation_volume(cell); + } + +double Mesh::cell_orientation_volume(const int* cell) const + { + // 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; } -std::vector Mesh::asShape(const std::vector& dx) const +void Mesh::compute_cell(int* cell, const double* coordinate) const { - std::vector temp; - for (int idx = 0; idx < shape_.size(); ++idx) + for (int dim = 0; dim < num_dim_; ++dim) { - temp.push_back(static_cast(std::ceil(dx[idx] / step_[idx]))); + cell[dim] = static_cast((coordinate[dim] - lower_[dim]) / step_[dim]); } - return temp; + + // shift position coordinate + assert(num_position_dim_ == 1); + cell[0] -= start_; } -std::vector Mesh::asLength(const std::vector& shape) const +void Mesh::length_as_shape(int* shape, const double* length) const { - std::vector temp; - for (int idx = 0; idx < shape_.size(); ++idx) + for (int dim = 0; dim < num_dim_; ++dim) { - temp.push_back(shape[idx] * step_[idx]); + shape[dim] = static_cast(std::ceil(length[dim] / step_[dim])); } - return temp; } -double Mesh::integrateSurface(const std::vector& idx, double j_lo, double j_hi) const +void Mesh::shape_as_length(double* length, const int* shape) const { - return area(idx) * j_lo - area(idx + 1) * j_hi; + for (int dim = 0; dim < num_dim_; ++dim) + { + length[dim] = shape[dim] * step_[dim]; + } } -double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const +int Mesh::num_position_coordinates() const { - return integrateSurface(idx, j(idx), j(idx + 1)); + return num_position_dim_; } -double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const +int Mesh::num_orientation_coordinates() const { - return integrateSurface(const std::vector& idx, j(idx), j(idx + 1)); + 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_; } -double Mesh::integrateVolume(std::vector idx, double f) const +BoundaryType Mesh::lower_boundary_condition() const { - return volume(idx) * f; + return lower_bc_; } -double Mesh::integrateVolume(const std::vector& idx, const DataView& f) const +BoundaryType Mesh::upper_boundary_condition() const { - return integrateVolume(const std::vector& idx, f(idx)); + return upper_bc_; } -double Mesh::integrateVolume(const std::vector& idx, const DataView& f) const +#if 0 +double Mesh::integrateSurface(const std::vector& idx, double j_lo, double j_hi) const { - return integrateVolume(const std::vector& idx, f(idx)); + return area(idx) * j_lo - area(idx + 1) * j_hi; + } + +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const + { + return integrateSurface(idx, j(idx), j(idx + 1)); + } + +double Mesh::integrateSurface(const std::vector& idx, const DataView& j) const + { + return integrateSurface(const std::vector& idx, j(idx), j(idx + 1)); } double Mesh::gradient(int idx, const DataView& f) const @@ -186,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 @@ -199,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/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 From 217824e26561632fbbda162f008f85a6b19a0985 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Wed, 9 Jul 2025 20:57:11 -0500 Subject: [PATCH 11/13] Continue to revise multidimensional arrays --- include/flyft/complex_field.h | 20 ++ include/flyft/data_array.h | 337 +++++++++++++++++++++++ include/flyft/data_layout.h | 199 ++++++++++++-- include/flyft/data_view.h | 490 ++++++++++++++++++++++------------ include/flyft/field.h | 200 +------------- include/flyft/mesh.h | 1 + include/flyft/parallel_mesh.h | 7 +- src/data_layout.cc | 208 +-------------- src/data_view.cc | 9 + src/parallel_mesh.cc | 315 +++++++++++++++------- 10 files changed, 1088 insertions(+), 698 deletions(-) create mode 100644 include/flyft/complex_field.h create mode 100644 include/flyft/data_array.h create mode 100644 src/data_view.cc 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..b741e12 --- /dev/null +++ b/include/flyft/data_array.h @@ -0,0 +1,337 @@ +#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 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 +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(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 f5c3ce2..c424c67 100644 --- a/include/flyft/data_layout.h +++ b/include/flyft/data_layout.h @@ -1,6 +1,7 @@ #ifndef FLYFT_DATA_LAYOUT_H_ #define FLYFT_DATA_LAYOUT_H_ +#include #include namespace flyft @@ -8,12 +9,17 @@ 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 + * 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 multi-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. + * 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: @@ -27,21 +33,6 @@ class DataLayout */ DataLayout(int num_dimensions, const int* shape); - //! Copy constructor - DataLayout(const DataLayout& other); - - //! Copy assignment operator - DataLayout& operator=(const DataLayout& other); - - //! Move constructor - DataLayout(const DataLayout&& other); - - //! Move assignment operator - DataLayout& operator=(const DataLayout&& other); - - //! Destructor - ~DataLayout(); - //! Map a multidimensional index to a one-dimensional index. /*! * \param multi_index Multidimensional index. @@ -66,7 +57,7 @@ class DataLayout */ void operator()(int* multi_index, size_t flat_index) const; - //! Map a one-dimensional index to a multidimensional index offset. + //! Map a one-dimensional index to a multidimensional index relative to an offset. /*! * \param multi_index Multidimensional index. * \param flat_index One-dimensional index. @@ -84,20 +75,176 @@ class DataLayout size_t size() const; //! Test if two layouts are equal. - bool operator==(const DataLayout& other) const; + template + bool operator==(const DataLayout& other) const; //! Test if two layouts are not equal. - bool operator!=(const DataLayout& other) const; + template + bool operator!=(const DataLayout& other) const; + + static int constexpr MAX_NUM_DIMENSIONS = N; //!< Maximum number of dimensions. private: int num_dimensions_; //!< Number of dimensions. - int* shape_; //!< Number of elements per dimension. + 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. - - //! Reset shape and dimensionality of layout - void reset(int num_dimensions, const int* shape); }; +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 4926c25..a08eb54 100644 --- a/include/flyft/data_view.h +++ b/include/flyft/data_view.h @@ -3,6 +3,8 @@ #include "flyft/data_layout.h" +#include +#include #include #include @@ -12,16 +14,19 @@ namespace flyft //! 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 multi-index, as well as a forward iterator. + * 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 multi-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 - * multi-index of 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 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 +template class DataView { public: @@ -29,227 +34,356 @@ class DataView using pointer = value_type*; //; //!< Layout type + class Iterator { public: - using value_type = DataView::value_type; //!< Data type - using pointer = DataView::pointer; //!< Pointer to data type - using reference = DataView::reference; //!< Reference to data type + 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 //! Empty constructor - Iterator() : Iterator(DataView()) {} + Iterator(); //! Construct from view. - explicit Iterator(const DataView& view) : Iterator(view, nullptr) {} - - //! Construct from view starting from multidimensional index. - Iterator(const DataView& view, int* current) : view_(view) - { - current_ = new int[view_.num_dimensions()]; - for (int dim = 0; dim < view_.num_dimensions(); ++dim) - { - current_[dim] = (current) ? current[dim] : 0; - } - } - - // non-copyable - Iterator(const Iterator& other) = delete; - Iterator& operator=(const Iterator& other) = delete; - - //! Move constructor - Iterator(const Iterator&& other) - : view_(std::move(other.view_)), current_(std::move(other.current_)) - { - } - - //! Move assignment operator. - Iterator& operator=(const Iterator&& other) - { - view_ = std::move(other.view_); - current_ = std::move(other.current_); - return *this; - } - - //! Destructor. - ~Iterator() - { - delete[] current_; - } - - reference operator*() const - { - return view_(current_); - } - - pointer operator->() const - { - return get(); - } - - pointer get() const - { - return &view_(current_); - } - - 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) - { - current_[dim] = 0; - ++current_[dim--]; - } - - return *this; - } - - bool operator==(const Iterator& other) const - { - return (get() == other.get()); - } - - bool operator!=(const Iterator& other) const - { - return !(*this == other); - } + /*! + * \param view Data view to iterate. + * + * This iterator starts at the zero multidimensional index. + */ + explicit Iterator(const DataView& view); + + //! 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); + + //! Dereference iterator. + reference operator*() const; + + //! Access current data as pointer. + pointer operator->() const; + + //! Get pointer to current data. + pointer get() const; + + //! Pre-increment. + /*! + * The last dimension of the multidimensional index is incremented by 1. + */ + Iterator& operator++(); + + //! Post-increment. + /*! + * The last dimension of the multidimensional index is incremented by 1. + */ + Iterator operator++(int); + + //! 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); + + //! 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); + + //! 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; + + //! 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 }; //! Empty constructor. - DataView() : data_(nullptr), layout_(DataLayout()), start_(nullptr), shape_(nullptr) {} + DataView(); //! Constructor. /*! * \param data Pointer to the data. * \param layout Layout of the data. */ - DataView(pointer data, const DataLayout& layout) : data_(data), layout_(layout) - { - start_ = new int[layout_.num_dimensions()]; - shape_ = new int[layout_.num_dimensions()]; - const auto layout_shape = layout_.shape(); - for (int dim = 0; dim < layout_.num_dimensions(); ++dim) - { - start_[dim] = 0; - shape_[dim] = layout_shape[dim]; - } - } + DataView(pointer data, const Layout& layout); //! 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. + * \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 DataLayout& layout, const int* start, const int* shape) - : data_(data), layout_(layout) - { - start_ = new int[layout_.num_dimensions()]; - shape_ = new int[layout_.num_dimensions()]; - for (int dim = 0; dim < layout_.num_dimensions(); ++dim) - { - start_[dim] = start[dim]; - shape_[dim] = shape[dim]; - } - } - - // non-copyable, non-movable - DataView(const DataView& other) = delete; - DataView& operator=(const DataView& other) = delete; - DataView(const DataView&& other) = delete; - DataView& operator=(const DataView&& other) = delete; - - //! Destructor. - ~DataView() - { - delete[] start_; - delete[] shape_; - } + DataView(pointer data, const Layout& layout, const int* start, const int* shape); //! 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 - { - return data_[layout_(multi_index, start_)]; - } + reference operator()(const int* multi_index) const; //! Number of dimensions. - int num_dimensions() const - { - return layout_.num_dimensions(); - } + int num_dimensions() const; //! Number of elements per dimension in the scope of the view. - const int* shape() const - { - return shape_; - } + const int* shape() const; //! Number of elements in the scope of the view. - size_t size() const - { - size_t size; - if (shape_) - { - size = 1; - for (int dim = 0; dim < layout_.num_dimensions(); ++dim) - { - size *= shape_[dim]; - } - } - else - { - size = 0; - } - return size; - } + size_t size() const; //! Test if view is not null /*! * \returns True if data pointer being viewed is not null. */ - explicit operator bool() const - { - return (data_ != nullptr); - } + explicit operator bool() const; - //! Start point of the data + //! Make an iterator to the start of the data being viewed. /*! * \return Iterator to start of data. */ - Iterator begin() const - { - return Iterator(*this); - } + Iterator begin() const; - //! End point of the data + //! Make and iterator to the end of the data being viewed. /*! * \return Iterator to end of data. */ - Iterator end() const - { - return Iterator(*this, shape_); - } + 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. private: - pointer data_; //!< Data. - DataLayout layout_; //!< Multidimensional array layout. - int* start_; //!< Multi-index for start of scoped view. - int* shape_; //!< Number of elements per dimension of scoped view. + 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) + { + start_[dim] = 0; + shape_[dim] = layout_shape[dim]; + } + } + +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) + { + assert(start_[dim] >= 0); + assert(shape_[dim] > 0); + assert(start_[dim] + shape_[dim] <= layout.shape()[dim]); + start_[dim] = start[dim]; + shape_[dim] = shape[dim]; + } + } + +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_) + { + size = 1; + for (int dim = 0; dim < layout_.num_dimensions(); ++dim) + { + size *= shape_[dim]; + } + } + else + { + 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); + } + +template +DataView::Iterator::Iterator(const DataView& view, const int* current) : view_(view) + { + assert(view); + for (int dim = 0; dim < view_.num_dimensions(); ++dim) + { + 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(); + } + +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) + { + current_[dim] = 0; + ++current_[dim--]; + } + + 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) + { + (*this)++; + } + return *this; + } + +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 #endif // FLYFT_DATA_VIEW_H_ diff --git a/include/flyft/field.h b/include/flyft/field.h index 9e77dbb..b24b841 100644 --- a/include/flyft/field.h +++ b/include/flyft/field.h @@ -1,203 +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(char num_dimensions, const int* shape) - : GenericField(num_dimensions, shape, nullptr) - { - } - - GenericField(char num_dimensions, const int* shape, const int* buffer_shape) - { - reshape(num_dimensions, shape, buffer_shape); - } - - // noncopyable / nonmovable - GenericField(const GenericField&) = delete; - GenericField(GenericField&&) = delete; - GenericField& operator=(const GenericField&) = delete; - GenericField& operator=(GenericField&&) = delete; - - ~GenericField() - { - delete[] data_; - delete[] shape_; - delete[] buffer_shape_; - delete[] full_shape_; - } - - T& operator()(const int* multi_index) - { - return view()(multi_index); - } - - const T& operator()(const int* multi_index) const - { - return const_view()(multi_index); - } - - const int* shape() const - { - return shape_; - } - - const int* buffer_shape() const - { - return buffer_shape_; - } - - const int* full_shape() const - { - return layout_.shape(); - } - - View view() - { - token_.stageAndCommit(); - return View(data_, layout_, buffer_shape_, shape_); - } - - ConstantView const_view() const - { - return ConstantView(static_cast(data_), layout_, buffer_shape_, shape_); - } - - View full_view() - { - token_.stageAndCommit(); - return View(data_, layout_); - } - - ConstantView const_full_view() const - { - return ConstantView(data_, layout_); - } - - void reshape(char num_dimensions, const int* shape, const int* buffer_shape) - { - if (num_dimensions != layout_.num_dimensions()) - { - delete[] shape_; - delete[] buffer_shape_; - delete[] full_shape_; - - shape_ = new int[num_dimensions] {}; - buffer_shape_ = new int[num_dimensions] {}; - full_shape_ = new int[num_dimensions] {}; - } - - bool same_shape = true; - bool same_buffer_shape = true; - for (char dim = 0; dim < num_dimensions; ++dim) - { - if (shape[dim] != shape_[dim]) - { - same_shape = false; - shape_[dim] = shape[dim]; - } - - if (buffer_shape[dim] != buffer_shape_[dim]) - { - same_buffer_shape = false; - buffer_shape_[dim] = buffer_shape[dim]; - } - - full_shape_[dim] = shape_[dim] + 2 * buffer_shape_[dim]; - } - - if (!(same_shape && same_buffer_shape)) - { - DataLayout layout(num_dimensions, full_shape_); - - // 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 = 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, 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_ && layout.size() > 0) - { - data_ = new T[layout.size()]; - std::fill(data_, data_ + layout.size(), T(0)); - } - layout_ = layout; - token_.stageAndCommit(); - } - } - - void setBuffer(const int* buffer_shape) - { - reshape(shape_, buffer_shape); - } - - void requestBuffer(const int* buffer_shape) - { - bool any_increase = false; - for (char dim = 0; !any_increase && dim < layout_.num_dimensions(); ++dim) - { - any_increase |= (buffer_shape[dim] > buffer_shape_[dim]); - } - - if (any_increase) - { - setBuffer(buffer_shape); - } - } - - private: - T* data_; - int* shape_; - int* buffer_shape_; - int* full_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 a2f1a9b..10a6879 100644 --- a/include/flyft/mesh.h +++ b/include/flyft/mesh.h @@ -52,6 +52,7 @@ class Mesh //! Move assignment operator. Mesh& operator=(const Mesh&& other); + //! Destructor virtual ~Mesh(); //! Slice mesh by position coordinate. diff --git a/include/flyft/parallel_mesh.h b/include/flyft/parallel_mesh.h index 0c6ecf5..a5555c6 100644 --- a/include/flyft/parallel_mesh.h +++ b/include/flyft/parallel_mesh.h @@ -28,7 +28,7 @@ class ParallelMesh 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 +38,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/src/data_layout.cc b/src/data_layout.cc index 56db9a9..f6ac194 100644 --- a/src/data_layout.cc +++ b/src/data_layout.cc @@ -1,214 +1,8 @@ #include "flyft/data_layout.h" -#include - namespace flyft { -DataLayout::DataLayout() : num_dimensions_(0), shape_(nullptr), size_(0) {} - -DataLayout::DataLayout(int num_dimensions, const int* shape) - { - reset(num_dimensions, shape); - } - -DataLayout::DataLayout(const DataLayout& other) - { - reset(other.num_dimensions_, other.shape_); - } - -DataLayout& DataLayout::operator=(const DataLayout& other) - { - if (this != &other) - { - reset(other.num_dimensions_, other.shape_); - } - return *this; - } - -DataLayout::DataLayout(const DataLayout&& other) - : num_dimensions_(std::move(other.num_dimensions_)), shape_(std::move(other.shape_)), - size_(std::move(other.size_)) - { - } - -DataLayout& DataLayout::operator=(const DataLayout&& other) - { - num_dimensions_ = std::move(other.num_dimensions_); - shape_ = std::move(other.shape_); - size_ = std::move(other.size_); - return *this; - } - -DataLayout::~DataLayout() - { - delete[] shape_; - } - -size_t DataLayout::operator()(const int* multi_index) const - { - size_t flat_index; - if (num_dimensions_ > 1) - { - // multiply out using row-major ordering - flat_index = 0; - size_t stride = 1; - for (int dim = num_dimensions_ - 1; dim >= 0; --dim) - { - flat_index += multi_index[dim] * stride; - stride *= shape_[dim]; - } - } - else - { - flat_index = multi_index[0]; - } - return flat_index; - } - -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 - size_t flat_index; - if (num_dimensions_ > 1) - { - // multiply out using row-major ordering - flat_index = 0; - size_t stride = 1; - for (int dim = num_dimensions_ - 1; dim >= 0; --dim) - { - flat_index += (offset[dim] + multi_index[dim]) * stride; - stride *= shape_[dim]; - } - } - else - { - flat_index = offset[0] + multi_index[0]; - } - return flat_index; - } - -void DataLayout::operator()(int* multi_index, size_t flat_index) const - { - if (num_dimensions_ > 1) - { - // repeatedly floor divide the flat index by stride associated with it - size_t idx = flat_index; - size_t stride = size_; - for (int dim = 0; dim < num_dimensions_ - 1; ++dim) - { - stride /= shape_[dim]; - multi_index[dim] = idx / stride; - idx -= multi_index[dim] * stride; - } - } - else - { - multi_index[0] = flat_index; - } - } - -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_ - 1; ++dim) - { - multi_index[dim] -= offset[dim]; - } - } - } - -int DataLayout::num_dimensions() const - { - return num_dimensions_; - } - -const int* DataLayout::shape() const - { - return shape_; - } - -size_t DataLayout::size() const - { - return size_; - } - -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) - { - // then, shape must match along all dimensions - if (shape_ && other.shape_) - { - for (int i = 0; i < num_dimensions_; ++i) - { - if (shape_[i] != other.shape_[i]) - { - same_layout = false; - } - } - } - // or, both pointers must be null - else if (!(!shape_ && !other.shape_)) - { - same_layout = false; - } - } - - return same_layout; - } - -bool DataLayout::operator!=(const DataLayout& other) const - { - return !(*this == other); - } - -void DataLayout::reset(int num_dimensions, const int* shape) - { - if (shape && num_dimensions > 0) - { - // free shape if it is allocated but the wrong size - if (shape_ && num_dimensions_ != num_dimensions) - { - delete[] shape_; - shape_ = nullptr; - } - - // (re-)allocate shape memory if needed - if (!shape_) - { - shape_ = new int[num_dimensions]; - } - - // always copy shape values and recompute size - num_dimensions_ = num_dimensions; - size_ = 1; - for (int i = 0; i < num_dimensions; ++i) - { - shape_[i] = shape[i]; - size_ *= shape[i]; - } - } - else - { - if (shape_) - { - delete[] shape_; - shape_ = nullptr; - } - num_dimensions_ = 0; - size_ = 0; - } - } +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/parallel_mesh.cc b/src/parallel_mesh.cc index d8f98e1..8445054 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); @@ -70,40 +72,43 @@ std::shared_ptr ParallelMesh::full() const 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 +140,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 +274,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 +303,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 +417,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 +438,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()); From feeace56b900629754576e6f96c07c522f553e97 Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 10 Jul 2025 14:10:29 -0500 Subject: [PATCH 12/13] Run pre-commit --- python/tests/conftest.py | 3 ++- python/tests/test_boublik_hard_sphere.py | 3 ++- python/tests/test_brownian_diffusive_flux.py | 3 ++- python/tests/test_composite_external_potential.py | 3 ++- python/tests/test_composite_flux.py | 3 ++- python/tests/test_composite_functional.py | 3 ++- python/tests/test_crank_nicolson_integrator.py | 3 ++- python/tests/test_explicit_euler_integrator.py | 3 ++- python/tests/test_exponential_wall_potential.py | 3 ++- python/tests/test_external_field.py | 3 ++- python/tests/test_field.py | 3 ++- python/tests/test_grand_potential.py | 3 ++- python/tests/test_hard_wall_potential.py | 3 ++- python/tests/test_harmonic_wall_potential.py | 3 ++- python/tests/test_implicit_euler_integrator.py | 3 ++- python/tests/test_lennard_jones_93_wall_potential.py | 3 ++- python/tests/test_mirror.py | 3 ++- python/tests/test_parameter.py | 3 ++- python/tests/test_picard_iteration.py | 3 ++- python/tests/test_rpy_diffusive_flux.py | 3 ++- python/tests/test_state.py | 3 ++- python/tests/test_virial_expansion.py | 3 ++- python/tests/test_white_bear.py | 3 ++- python/tests/test_white_bear_mark_ii.py | 3 ++- python/validate/binary_mixture.py | 3 ++- python/validate/brownian_diffusion.py | 3 ++- python/validate/hard_sphere_evaporation.py | 3 ++- python/validate/hard_sphere_slit.py | 3 ++- python/validate/stratification.py | 3 ++- 29 files changed, 58 insertions(+), 29 deletions(-) diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 4adb701..04d65ed 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -1,7 +1,8 @@ -import flyft import pytest from pytest_lazy_fixtures import lf as lazy_fixture +import flyft + @pytest.fixture def ig(): diff --git a/python/tests/test_boublik_hard_sphere.py b/python/tests/test_boublik_hard_sphere.py index 98c4463..0b7c3cb 100644 --- a/python/tests/test_boublik_hard_sphere.py +++ b/python/tests/test_boublik_hard_sphere.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + def fex_cs(eta, d): """Carnahan-Starling free-energy density of hard spheres""" diff --git a/python/tests/test_brownian_diffusive_flux.py b/python/tests/test_brownian_diffusive_flux.py index e55ad57..f7099c1 100644 --- a/python/tests/test_brownian_diffusive_flux.py +++ b/python/tests/test_brownian_diffusive_flux.py @@ -1,8 +1,9 @@ -import flyft import numpy as np import pytest from pytest_lazy_fixtures import lf as lazy_fixture +import flyft + @pytest.fixture def cartesian_mesh_grand(): diff --git a/python/tests/test_composite_external_potential.py b/python/tests/test_composite_external_potential.py index e075944..9515921 100644 --- a/python/tests/test_composite_external_potential.py +++ b/python/tests/test_composite_external_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def comp(): diff --git a/python/tests/test_composite_flux.py b/python/tests/test_composite_flux.py index b92b7d0..f1785eb 100644 --- a/python/tests/test_composite_flux.py +++ b/python/tests/test_composite_flux.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def comp(): diff --git a/python/tests/test_composite_functional.py b/python/tests/test_composite_functional.py index ec070a1..ebd37f6 100644 --- a/python/tests/test_composite_functional.py +++ b/python/tests/test_composite_functional.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def comp(): diff --git a/python/tests/test_crank_nicolson_integrator.py b/python/tests/test_crank_nicolson_integrator.py index 2a9153f..dfcba14 100644 --- a/python/tests/test_crank_nicolson_integrator.py +++ b/python/tests/test_crank_nicolson_integrator.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def cn(): diff --git a/python/tests/test_explicit_euler_integrator.py b/python/tests/test_explicit_euler_integrator.py index 1e0e3e6..2889542 100644 --- a/python/tests/test_explicit_euler_integrator.py +++ b/python/tests/test_explicit_euler_integrator.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def euler(): diff --git a/python/tests/test_exponential_wall_potential.py b/python/tests/test_exponential_wall_potential.py index 18573b1..5da332b 100644 --- a/python/tests/test_exponential_wall_potential.py +++ b/python/tests/test_exponential_wall_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def ew(): diff --git a/python/tests/test_external_field.py b/python/tests/test_external_field.py index f896a48..b07c5d0 100644 --- a/python/tests/test_external_field.py +++ b/python/tests/test_external_field.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def spherical_mesh_grand(): diff --git a/python/tests/test_field.py b/python/tests/test_field.py index c7f3d78..2409aff 100644 --- a/python/tests/test_field.py +++ b/python/tests/test_field.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def field(): diff --git a/python/tests/test_grand_potential.py b/python/tests/test_grand_potential.py index bbf9897..768cd2c 100644 --- a/python/tests/test_grand_potential.py +++ b/python/tests/test_grand_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + from .test_ideal_gas import f_ig, mu_ig from .test_rosenfeld_fmt import fex_py, muex_py diff --git a/python/tests/test_hard_wall_potential.py b/python/tests/test_hard_wall_potential.py index 5730945..697d82e 100644 --- a/python/tests/test_hard_wall_potential.py +++ b/python/tests/test_hard_wall_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def hw(): diff --git a/python/tests/test_harmonic_wall_potential.py b/python/tests/test_harmonic_wall_potential.py index a700c2a..7d6dbd2 100644 --- a/python/tests/test_harmonic_wall_potential.py +++ b/python/tests/test_harmonic_wall_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def hw(): diff --git a/python/tests/test_implicit_euler_integrator.py b/python/tests/test_implicit_euler_integrator.py index dc7bc7e..aa39865 100644 --- a/python/tests/test_implicit_euler_integrator.py +++ b/python/tests/test_implicit_euler_integrator.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def euler(): diff --git a/python/tests/test_lennard_jones_93_wall_potential.py b/python/tests/test_lennard_jones_93_wall_potential.py index ca70f83..efa5692 100644 --- a/python/tests/test_lennard_jones_93_wall_potential.py +++ b/python/tests/test_lennard_jones_93_wall_potential.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def lj(): diff --git a/python/tests/test_mirror.py b/python/tests/test_mirror.py index 1828282..f721114 100644 --- a/python/tests/test_mirror.py +++ b/python/tests/test_mirror.py @@ -1,6 +1,7 @@ -import flyft.mirror import pytest +import flyft.mirror + class _A: def __init__(self, x, y): diff --git a/python/tests/test_parameter.py b/python/tests/test_parameter.py index 3e91ea5..82631cc 100644 --- a/python/tests/test_parameter.py +++ b/python/tests/test_parameter.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + class TimeParameter(flyft.CustomParameter): def __call__(self, state): diff --git a/python/tests/test_picard_iteration.py b/python/tests/test_picard_iteration.py index b5ad3de..5bc1cf5 100644 --- a/python/tests/test_picard_iteration.py +++ b/python/tests/test_picard_iteration.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + from .test_ideal_gas import mu_ig from .test_rosenfeld_fmt import muex_py diff --git a/python/tests/test_rpy_diffusive_flux.py b/python/tests/test_rpy_diffusive_flux.py index 89e433c..71d4524 100644 --- a/python/tests/test_rpy_diffusive_flux.py +++ b/python/tests/test_rpy_diffusive_flux.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + @pytest.fixture def spherical_mesh_grand(): diff --git a/python/tests/test_state.py b/python/tests/test_state.py index a1aed37..24924a6 100644 --- a/python/tests/test_state.py +++ b/python/tests/test_state.py @@ -1,6 +1,7 @@ -import flyft import pytest +import flyft + def test_init(state): assert isinstance(state.mesh, flyft.state.ParallelMesh) diff --git a/python/tests/test_virial_expansion.py b/python/tests/test_virial_expansion.py index 3464b95..4b5d0e1 100644 --- a/python/tests/test_virial_expansion.py +++ b/python/tests/test_virial_expansion.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + def f_ex(B, rho): """Excess free-energy density of virial expansion""" diff --git a/python/tests/test_white_bear.py b/python/tests/test_white_bear.py index 6ef571b..2924ae2 100644 --- a/python/tests/test_white_bear.py +++ b/python/tests/test_white_bear.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + def fex_cs(eta, v): """Carnahan--Starling free-energy density of hard spheres""" diff --git a/python/tests/test_white_bear_mark_ii.py b/python/tests/test_white_bear_mark_ii.py index 48ee0ec..73feed1 100644 --- a/python/tests/test_white_bear_mark_ii.py +++ b/python/tests/test_white_bear_mark_ii.py @@ -1,7 +1,8 @@ -import flyft import numpy as np import pytest +import flyft + def fex_cs(eta, v): """Carnahan--Starling free-energy density of hard spheres""" diff --git a/python/validate/binary_mixture.py b/python/validate/binary_mixture.py index 28e6620..4cdae6b 100644 --- a/python/validate/binary_mixture.py +++ b/python/validate/binary_mixture.py @@ -1,6 +1,7 @@ -import flyft import numpy as np +import flyft + L = 50.0 dx = 0.2 diameters = {"A": 2.0, "B": 4.0} diff --git a/python/validate/brownian_diffusion.py b/python/validate/brownian_diffusion.py index 0fc5d5a..8183677 100644 --- a/python/validate/brownian_diffusion.py +++ b/python/validate/brownian_diffusion.py @@ -1,6 +1,7 @@ -import flyft import numpy as np +import flyft + L = 10.0 dx = 0.05 diameters = {"A": 1.0} diff --git a/python/validate/hard_sphere_evaporation.py b/python/validate/hard_sphere_evaporation.py index be85854..af01e9e 100644 --- a/python/validate/hard_sphere_evaporation.py +++ b/python/validate/hard_sphere_evaporation.py @@ -1,6 +1,7 @@ -import flyft import numpy as np +import flyft + L = 10.0 dx = 0.05 diameters = {"A": 1.0} diff --git a/python/validate/hard_sphere_slit.py b/python/validate/hard_sphere_slit.py index bb2a014..73d6fdd 100644 --- a/python/validate/hard_sphere_slit.py +++ b/python/validate/hard_sphere_slit.py @@ -1,6 +1,7 @@ -import flyft import numpy as np +import flyft + L = 16.0 diameters = {"A": 1.0} etas = {"A": 0.4257} diff --git a/python/validate/stratification.py b/python/validate/stratification.py index ee04a00..a8cdb6f 100644 --- a/python/validate/stratification.py +++ b/python/validate/stratification.py @@ -1,6 +1,7 @@ -import flyft import numpy as np +import flyft + L = 38.5 dx = 0.05 diameters = {"A": 1.0, "B": 4.0} From 7569334bf3b670cc15c867386f73f720dd57f16e Mon Sep 17 00:00:00 2001 From: Michael Howard Date: Thu, 10 Jul 2025 20:57:03 -0500 Subject: [PATCH 13/13] Update state --- include/flyft/data_array.h | 11 ++++- include/flyft/parallel_mesh.h | 9 ++++ include/flyft/state.h | 86 +++++++++++++++++++++++++++++----- src/CMakeLists.txt | 59 +++++++++++------------ src/parallel_mesh.cc | 15 ++++++ src/state.cc | 88 ++++++++++++----------------------- 6 files changed, 169 insertions(+), 99 deletions(-) diff --git a/include/flyft/data_array.h b/include/flyft/data_array.h index b741e12..d671e86 100644 --- a/include/flyft/data_array.h +++ b/include/flyft/data_array.h @@ -86,6 +86,9 @@ class DataArray : public TrackedObject */ const T& operator()(const int* multi_index) const; + //! Number of dimensions. + int num_dimensions() const; + //! Number of elements per dimension. const int* shape() const; @@ -189,6 +192,12 @@ const T& DataArray::operator()(const int* multi_index) const return const_view()(multi_index); } +template +int DataArray::num_dimensions() const + { + return layout_.num_dimensions(); + } + template const int* DataArray::shape() const { @@ -307,7 +316,7 @@ void DataArray::reshape(int num_dimensions, const int* shape, const int* b template void DataArray::setBuffer(const int* buffer_shape) { - reshape(shape_, buffer_shape); + reshape(layout_.num_dimensions(), shape_, buffer_shape); } template diff --git a/include/flyft/parallel_mesh.h b/include/flyft/parallel_mesh.h index a5555c6..3584a47 100644 --- a/include/flyft/parallel_mesh.h +++ b/include/flyft/parallel_mesh.h @@ -25,6 +25,15 @@ class ParallelMesh std::shared_ptr 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; 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/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/parallel_mesh.cc b/src/parallel_mesh.cc index 8445054..a2152d3 100644 --- a/src/parallel_mesh.cc +++ b/src/parallel_mesh.cc @@ -70,6 +70,21 @@ 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 comm_->size(); 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()); } } }