From b406f5279ab5172736f43898da291a7a3aa3419c Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 18 Mar 2021 14:18:41 -0700 Subject: [PATCH 01/43] MultiFab: Basics --- src/Base/CMakeLists.txt | 1 + src/Base/MultiFab.cpp | 171 ++++++++++++++++++++++++++++++++++++++++ src/pyAMReX.cpp | 2 + 3 files changed, 174 insertions(+) create mode 100644 src/Base/MultiFab.cpp diff --git a/src/Base/CMakeLists.txt b/src/Base/CMakeLists.txt index 024f8f0d..a7297c8f 100644 --- a/src/Base/CMakeLists.txt +++ b/src/Base/CMakeLists.txt @@ -4,4 +4,5 @@ target_sources(pyAMReX Box.cpp Dim3.cpp IntVect.cpp + MultiFab.cpp ) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp new file mode 100644 index 00000000..966f6643 --- /dev/null +++ b/src/Base/MultiFab.cpp @@ -0,0 +1,171 @@ +/* Copyright 2021 The AMReX Community + * + * Authors: Axel Huebl, ... + * License: BSD-3-Clause-LBNL + */ +#include +#include + +#include +#include + +#include + +namespace py = pybind11; +using namespace amrex; + + +void init_MultiFab(py::module &m) { + py::class_< MultiFab >(m, "MultiFab") + .def("__repr__", + [](MultiFab const & mf) { + return ""; + } + ) + + /* Constructors */ + .def(py::init< >()) + .def(py::init< const BoxArray&, const DistributionMapping&, int, int, + MFInfo const &, FabFactoryconst & >()) + .def(py::init< const BoxArray&, const DistributionMapping&, int, + IntVect const&, + MFInfo const&, FabFactory const & >()) + .def(py::init< MultiFab const&, MakeType, int, int >()) + + /* delayed defines */ + .def("define", + py::overload_cast< const BoxArray&, const DistributionMapping&, int, int, + MFInfo const &, FabFactory const & + >(&MultiFab::define)) + .def("define", + py::overload_cast< const BoxArray&, const DistributionMapping&, int, + IntVect const&, MFInfo const &, FabFactory const & + >(&MultiFab::define)) + + /* sizes, etc. */ + .def("min", + py::overload_cast< int, int, bool >(&MultiFab::min, py::const_)) + .def("min", + py::overload_cast< Box const &, int, int, bool >(&MultiFab::min, py::const_)) + .def("max", + py::overload_cast< int, int, bool >(&MultiFab::max, py::const_)) + .def("max", + py::overload_cast< Box const &, int, int, bool >(&MultiFab::max, py::const_)) + .def("minIndex", &MultiFab::minIndex) + .def("maxIndex", &MultiFab::maxIndex) + + /* norms */ + .def("norm0", py::overload_cast< int, int, bool, bool >(&MultiFab::norm0, py::const_)) + //.def("norm0", py::overload_cast< iMultiFab const &, int, int, bool >(&MultiFab::norm0, py::const_)) + + .def("norminf", py::overload_cast< int, int, bool, bool >(&MultiFab::norminf, py::const_)) + //.def("norminf", py::overload_cast< iMultiFab const &, int, int, bool >(&MultiFab::norminf, py::const_)) + + .def("norm1", py::overload_cast< int, Periodicity const&, bool >(&MultiFab::norm1, py::const_)) + .def("norm1", py::overload_cast< int, int, bool >(&MultiFab::norm1, py::const_)) + .def("norm1", py::overload_cast< Vector const &, int, bool >(&MultiFab::norm1, py::const_)) + + .def("norm2", py::overload_cast< int >(&MultiFab::norm2, py::const_)) + .def("norm2", py::overload_cast< int, Periodicity const& >(&MultiFab::norm2, py::const_)) + .def("norm2", py::overload_cast< Vector const & >(&MultiFab::norm2, py::const_)) + + /* simple math */ + .def("sum", &MultiFab::sum) + + .def("plus", py::overload_cast< Real, int >(&MultiFab::plus)) + .def("plus", py::overload_cast< Real, int, int, int >(&MultiFab::plus)) + .def("plus", py::overload_cast< Real, Box const &, int >(&MultiFab::plus)) + .def("plus", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::plus)) + .def("plus", py::overload_cast< MultiFab const &, int, int, int >(&MultiFab::plus)) + + .def("minus", py::overload_cast< MultiFab const &, int, int, int >(&MultiFab::minus)) + + // renamed: ImportError: overloading a method with both static and instance methods is not supported + .def("divi", py::overload_cast< MultiFab const &, int, int, int >(&MultiFab::divide)) + + .def("mult", py::overload_cast< Real, int >(&MultiFab::mult)) + .def("mult", py::overload_cast< Real, int, int, int >(&MultiFab::mult)) + .def("mult", py::overload_cast< Real, Box const &, int >(&MultiFab::mult)) + .def("mult", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::mult)) + + .def("invert", py::overload_cast< Real, int >(&MultiFab::invert)) + .def("invert", py::overload_cast< Real, int, int, int >(&MultiFab::invert)) + .def("invert", py::overload_cast< Real, Box const &, int >(&MultiFab::invert)) + .def("invert", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::invert)) + + .def("negate", py::overload_cast< int >(&MultiFab::negate)) + .def("negate", py::overload_cast< int, int, int >(&MultiFab::negate)) + .def("negate", py::overload_cast< Box const &, int >(&MultiFab::negate)) + .def("negate", py::overload_cast< Box const &, int, int, int >(&MultiFab::negate)) + + .def("sum_boundary", py::overload_cast< Periodicity const & >(&MultiFab::SumBoundary)) + .def("sum_boundary", py::overload_cast< int, int, Periodicity const & >(&MultiFab::SumBoundary)) + .def("sum_boundary", py::overload_cast< int, int, IntVect const&, Periodicity const & >(&MultiFab::SumBoundary)) + + /* static (standalone) simple math functions */ + .def_static("dot", py::overload_cast< MultiFab const &, int, MultiFab const &, int, int, int, bool >(&MultiFab::Dot)) + .def_static("dot", py::overload_cast< MultiFab const &, int, int, int, bool >(&MultiFab::Dot)) + //.def_static("dot", py::overload_cast< iMultiFab const&, const MultiFab&, int, MultiFab const&, int, int, int, bool >(&MultiFab::Dot)) + + .def_static("add", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, int >(&MultiFab::Add)) + .def_static("add", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Add)) + + .def_static("subtract", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, int >(&MultiFab::Subtract)) + .def_static("subtract", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Subtract)) + + .def_static("multiply", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, int >(&MultiFab::Multiply)) + .def_static("multiply", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Multiply)) + + .def_static("divide", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, int >(&MultiFab::Divide)) + .def_static("divide", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Divide)) + + .def_static("copy", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, int >(&MultiFab::Copy)) + .def_static("copy", py::overload_cast< MultiFab &, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Copy)) + + .def_static("swap", py::overload_cast< MultiFab &, MultiFab &, int, int, int, int >(&MultiFab::Swap)) + .def_static("swap", py::overload_cast< MultiFab &, MultiFab &, int, int, int, IntVect const & >(&MultiFab::Swap)) + + .def_static("saxpy", py::overload_cast< MultiFab &, Real, MultiFab const &, int, int, int, int >(&MultiFab::Saxpy)) + .def_static("saxpy", py::overload_cast< MultiFab &, Real, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Saxpy)) + + .def_static("xpay", py::overload_cast< MultiFab &, Real, MultiFab const &, int, int, int, int >(&MultiFab::Xpay)) + .def_static("xpay", py::overload_cast< MultiFab &, Real, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::Xpay)) + + .def_static("lin_comb", py::overload_cast< MultiFab &, Real, MultiFab const &, int, Real, MultiFab const &, int, int, int, int >(&MultiFab::LinComb)) + .def_static("lin_comb", py::overload_cast< MultiFab &, Real, MultiFab const &, int, Real, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::LinComb)) + + .def_static("add_product", py::overload_cast< MultiFab &, MultiFab const &, int, MultiFab const &, int, int, int, int >(&MultiFab::AddProduct)) + .def_static("add_product", py::overload_cast< MultiFab &, MultiFab const &, int, MultiFab const &, int, int, int, IntVect const & >(&MultiFab::AddProduct)) + + /* simple data validity checks */ + .def("contains_nan", py::overload_cast< bool >(&MultiFab::contains_nan, py::const_)) + .def("contains_nan", py::overload_cast< int, int, int, bool >(&MultiFab::contains_nan, py::const_)) + .def("contains_nan", py::overload_cast< int, int, IntVect const &, bool >(&MultiFab::contains_nan, py::const_)) + + .def("contains_inf", py::overload_cast< bool >(&MultiFab::contains_inf, py::const_)) + .def("contains_inf", py::overload_cast< int, int, int, bool >(&MultiFab::contains_inf, py::const_)) + .def("contains_inf", py::overload_cast< int, int, IntVect const &, bool >(&MultiFab::contains_inf, py::const_)) + + /* masks & ownership */ + // TODO: + // - OverlapMask -> std::unique_ptr + // - OwnerMask -> std::unique_ptr + + /* Syncs */ + .def("average_sync", &MultiFab::AverageSync) + .def("weighted_sync", &MultiFab::WeightedSync) + .def("override_sync", py::overload_cast< Periodicity const & >(&MultiFab::OverrideSync)) + //.def("override_sync", py::overload_cast< iMultiFab const &, Periodicity const & >(&MultiFab::OverrideSync)) + + /* Init & Finalize */ + .def_static("initialize", &MultiFab::Initialize ) + .def_static("finalize", &MultiFab::Finalize ) + + /* data access in Box index space */ + //.def("__iter__") + //.def("array", [](MultiFab & mf, MFIter mfi) { + // // return Array4Real... + //}) + ; +} diff --git a/src/pyAMReX.cpp b/src/pyAMReX.cpp index f222863e..32ee9ae3 100644 --- a/src/pyAMReX.cpp +++ b/src/pyAMReX.cpp @@ -20,6 +20,7 @@ void init_AMReX(py::module&); void init_Box(py::module &); void init_Dim3(py::module&); void init_IntVect(py::module &); +void init_MultiFab(py::module &); PYBIND11_MODULE(amrex_pybind, m) { m.doc() = R"pbdoc( @@ -38,6 +39,7 @@ PYBIND11_MODULE(amrex_pybind, m) { init_Dim3(m); init_IntVect(m); init_Box(m); + init_MultiFab(m); // API runtime version // note PEP-440 syntax: x.y.zaN but x.y.z.devN From 1f3c4c200268cab7609a0c92d4269269b2428d5d Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 6 Oct 2021 11:33:38 -0700 Subject: [PATCH 02/43] MultiFab.__iter__ --- src/Base/MultiFab.cpp | 38 ++++++++++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 966f6643..4dc276ae 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -9,11 +9,37 @@ #include #include +#include #include namespace py = pybind11; using namespace amrex; +namespace { + /** STL-like iteration over amrex::MFIter + * + * The amrex::MFIter interface is currently a bit too tricky to implement + * std::begin() and std::end() safely with OpenMP threading. + */ + class MFIterWrapper { + std::shared_ptr mfi; + public: + explicit MFIterWrapper(const MultiFab& mf) { + // For tiling support (OpenMP/thread pools) later on: + // MFIter mfi(mf, TilingIfNotGPU()); + mfi = std::make_shared(mf); + } + std::shared_ptr operator*() { return mfi; } + std::shared_ptr operator*() const { return mfi; } + MFIterWrapper& operator++() { ++(*mfi); return *this; } + }; + + class ValidSentinel {}; + + bool operator==(MFIterWrapper const& it, ValidSentinel const&) { + return (*it)->isValid(); + } +} void init_MultiFab(py::module &m) { py::class_< MultiFab >(m, "MultiFab") @@ -163,9 +189,13 @@ void init_MultiFab(py::module &m) { .def_static("finalize", &MultiFab::Finalize ) /* data access in Box index space */ - //.def("__iter__") - //.def("array", [](MultiFab & mf, MFIter mfi) { - // // return Array4Real... - //}) + .def("__iter__", + [](const MultiFab& mf) { + return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); + }, + /* Essential: keep object alive while iterator exists */ + py::keep_alive<0, 1>() + //py::return_value_policy::reference_internal + ) ; } From c990e0887aaf286245901e052d4c612ea188a182 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 7 Oct 2021 11:36:04 -0700 Subject: [PATCH 03/43] More MultiFab Updates --- src/Base/MultiFab.cpp | 56 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 53 insertions(+), 3 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 4dc276ae..b6ce8a03 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -7,6 +7,10 @@ #include #include +#include +#include +#include +#include #include #include @@ -42,7 +46,11 @@ namespace { } void init_MultiFab(py::module &m) { - py::class_< MultiFab >(m, "MultiFab") + py::class_< FabArrayBase >(m, "FabArrayBase"); + py::class_< FArrayBox >(m, "FArrayBox"); + py::class_< FabArray, FabArrayBase >(m, "FabArray_FArrayBox"); + + py::class_< MultiFab /*, FabArray*/ >(m, "MultiFab") .def("__repr__", [](MultiFab const & mf) { return "()) + //.def(py::init< MultiFab && >()) + .def(py::init< const BoxArray&, const DistributionMapping&, int, int >()) + .def(py::init< const BoxArray&, const DistributionMapping&, int, int, + MFInfo const & >()) .def(py::init< const BoxArray&, const DistributionMapping&, int, int, - MFInfo const &, FabFactoryconst & >()) + MFInfo const &, FabFactory const & >()) + + .def(py::init< const BoxArray&, const DistributionMapping&, int, + IntVect const& >()) + .def(py::init< const BoxArray&, const DistributionMapping&, int, + IntVect const&, + MFInfo const& >()) .def(py::init< const BoxArray&, const DistributionMapping&, int, IntVect const&, MFInfo const&, FabFactory const & >()) + .def(py::init< MultiFab const&, MakeType, int, int >()) /* delayed defines */ @@ -69,15 +88,42 @@ void init_MultiFab(py::module &m) { IntVect const&, MFInfo const &, FabFactory const & >(&MultiFab::define)) + /* setters */ + .def("set_val", + py::overload_cast< Real >(&MultiFab::setVal)) + .def("set_val", + [](MultiFab & mf, Real val, int comp, int num_comp) { + mf.setVal(val, comp, num_comp); + } + ) + .def("set_val", + py::overload_cast< Real, int, int, int >(&MultiFab::setVal)) + .def("set_val", + py::overload_cast< Real, int, int, IntVect const & >(&MultiFab::setVal)) + + .def("is_cell_centered", &MultiFab::is_cell_centered) + .def("is_nodal", + py::overload_cast< >(&MultiFab::is_nodal, py::const_)) + .def("is_nodal", + py::overload_cast< int >(&MultiFab::is_nodal, py::const_)) + + .def_property_readonly("num_comp", &MultiFab::nComp) + /* sizes, etc. */ + .def("min", + [](MultiFab const & mf, int comp) { mf.min(comp); }) .def("min", py::overload_cast< int, int, bool >(&MultiFab::min, py::const_)) .def("min", py::overload_cast< Box const &, int, int, bool >(&MultiFab::min, py::const_)) + + .def("max", + [](MultiFab const & mf, int comp) { mf.max(comp); }) .def("max", py::overload_cast< int, int, bool >(&MultiFab::max, py::const_)) .def("max", py::overload_cast< Box const &, int, int, bool >(&MultiFab::max, py::const_)) + .def("minIndex", &MultiFab::minIndex) .def("maxIndex", &MultiFab::maxIndex) @@ -99,6 +145,10 @@ void init_MultiFab(py::module &m) { /* simple math */ .def("sum", &MultiFab::sum) + .def("abs", + [](MultiFab & mf, int comp, int num_comp) { mf.abs(comp, num_comp); }) + .def("abs", py::overload_cast< int, int, int >(&MultiFab::abs)) + .def("plus", py::overload_cast< Real, int >(&MultiFab::plus)) .def("plus", py::overload_cast< Real, int, int, int >(&MultiFab::plus)) .def("plus", py::overload_cast< Real, Box const &, int >(&MultiFab::plus)) @@ -193,7 +243,7 @@ void init_MultiFab(py::module &m) { [](const MultiFab& mf) { return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); }, - /* Essential: keep object alive while iterator exists */ + // Essential: keep object alive while iterator exists py::keep_alive<0, 1>() //py::return_value_policy::reference_internal ) From bd50c51d3ac80e1b3d37e5b130d73bbbd3cb357d Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 7 Oct 2021 11:39:16 -0700 Subject: [PATCH 04/43] Add: BoxArray, DistributionMapping --- src/Base/BoxArray.cpp | 269 +++++++++++++++++++++++++++++++ src/Base/CMakeLists.txt | 2 + src/Base/DistributionMapping.cpp | 79 +++++++++ src/pyAMReX.cpp | 7 + 4 files changed, 357 insertions(+) create mode 100644 src/Base/BoxArray.cpp create mode 100644 src/Base/DistributionMapping.cpp diff --git a/src/Base/BoxArray.cpp b/src/Base/BoxArray.cpp new file mode 100644 index 00000000..e63ab6f6 --- /dev/null +++ b/src/Base/BoxArray.cpp @@ -0,0 +1,269 @@ +/* Copyright 2021 The AMReX Community + * + * Authors: Axel Huebl + * License: BSD-3-Clause-LBNL + */ +#include +#include + +#include +#include +#include + +#include + +namespace py = pybind11; +using namespace amrex; + + +void init_BoxArray(py::module &m) { + /* A collection of Boxes stored in an Array. It is a + * reference-counted concrete class, not a polymorphic one; i.e. you + * cannot use any of the List member functions with a BoxList + */ + py::class_< BoxArray >(m, "BoxArray") + .def("__repr__", + [](BoxArray const & ba) { + std::stringstream s; + s << ba.size(); + return ""; + } + ) + + //! Construct an empty BoxArray + .def(py::init<>()) + //.def(py::init< BoxArray const& >()) + //.def(py::init< BoxArray const&& >()) + + //! Construct a BoxArray from an array of Boxes of size nbox. + .def(py::init< Box const & >()) + //! Construct a BoxArray from an array of Boxes of size nbox. + .def(py::init< const Box*, int >()) + + /* + //! Construct a BoxArray from a BoxList. + explicit BoxArray (const BoxList& bl); + explicit BoxArray (BoxList&& bl) noexcept; + + BoxArray (const BoxArray& rhs, const BATransformer& trans); + + BoxArray (BoxList&& bl, IntVect const& max_grid_size); + */ + + .def_property_readonly("size", &BoxArray::size) + .def_property_readonly("capacity", &BoxArray::capacity) + .def_property_readonly("empty", &BoxArray::empty) + .def_property_readonly("numPts", &BoxArray::numPts) + .def_property_readonly("d_numPts", &BoxArray::d_numPts) +/* + .def_property("type", + py::overload_cast<>(&BoxArray::type, py::const_), + &Box::setType) + + .def_property_readonly("length", + py::overload_cast<>(&Box::length, py::const_)) + .def_property_readonly("is_empty", &Box::isEmpty) +*/ + + .def("define", + py::overload_cast< Box const & >(&BoxArray::define)) + //.def("define", + // py::overload_cast< BoxList const & >(&BoxArray::define)) + //.def("define", + // py::overload_cast< BoxList&& >(&BoxArray::define)) + + .def("clear", &BoxArray::clear) + .def("resize", &BoxArray::resize) + + .def("cell_equal", &BoxArray::CellEqual) + + .def("max_size", + py::overload_cast< int >(&BoxArray::maxSize)) + .def("max_size", + py::overload_cast< IntVect const& >(&BoxArray::maxSize)) + + .def("refine", + py::overload_cast< int >(&BoxArray::refine)) + .def("refine", + py::overload_cast< IntVect const & >(&BoxArray::refine)) + + //! Coarsen each Box in the BoxArray to the specified ratio. + .def("coarsen", + py::overload_cast< IntVect const & >(&BoxArray::coarsen)) + .def("coarsen", + py::overload_cast< int >(&BoxArray::coarsen)) + + //! Coarsen each Box in the BoxArray to the specified ratio. + .def("coarsenable", + py::overload_cast< int, int >(&BoxArray::coarsenable, py::const_)) + .def("coarsenable", + py::overload_cast< IntVect const &, int >(&BoxArray::coarsenable, py::const_)) + .def("coarsenable", + py::overload_cast< IntVect const &, IntVect const & >(&BoxArray::coarsenable, py::const_)) + +/* + //! Grow and then coarsen each Box in the BoxArray. + BoxArray& growcoarsen (int n, const IntVect& refinement_ratio); + BoxArray& growcoarsen (IntVect const& ngrow, const IntVect& refinement_ratio); + + //! Grow each Box in the BoxArray by the specified amount. + BoxArray& grow (int n); + + //! Grow each Box in the BoxArray by the specified amount. + BoxArray& grow (const IntVect& iv); + + //! \brief Grow each Box in the BoxArray on the low and high ends + //! by n_cell cells in the idir direction. + BoxArray& grow (int idir, int n_cell); + + //! \brief Grow each Box in the BoxArray on the low end + /?! by n_cell cells in the idir direction. + BoxArray& growLo (int idir, int n_cell); + + //! \brief Grow each Box in the BoxArray on the high end + //! by n_cell cells in the idir direction. + BoxArray& growHi (int idir, int n_cell); + //! \brief Apply surroundingNodes(Box) to each Box in BoxArray. + //! See the documentation of Box for details. + BoxArray& surroundingNodes (); + //! + //! \brief Apply surroundingNodes(Box,int) to each Box in + //! BoxArray. See the documentation of Box for details. + BoxArray& surroundingNodes (int dir); + + //! Apply Box::enclosedCells() to each Box in the BoxArray. + BoxArray& enclosedCells (); + + //! Apply Box::enclosedCells(int) to each Box in the BoxArray. + BoxArray& enclosedCells (int dir); + + //! Apply Box::convert(IndexType) to each Box in the BoxArray. + BoxArray& convert (IndexType typ); + + BoxArray& convert (const IntVect& typ); + + //! Apply function (*fp)(Box) to each Box in the BoxArray. + BoxArray& convert (Box (*fp)(const Box&)); + + //! Apply Box::shift(int,int) to each Box in the BoxArray. + BoxArray& shift (int dir, int nzones); + + //! Apply Box::shift(const IntVect &iv) to each Box in the BoxArray. + BoxArray& shift (const IntVect &iv); + + //! Set element i in this BoxArray to Box ibox. + void set (int i, const Box& ibox); + + //! Return element index of this BoxArray. + Box operator[] (int index) const noexcept { + return m_bat(m_ref->m_abox[index]); + } + + //! Return element index of this BoxArray. + Box operator[] (const MFIter& mfi) const noexcept; + + //! Return element index of this BoxArray. + Box get (int index) const noexcept { return operator[](index); } + + //! Return cell-centered box at element index of this BoxArray. + Box getCellCenteredBox (int index) const noexcept { + return m_bat.coarsen(m_ref->m_abox[index]); + } + + //! \brief Return true if Box is valid and they all have the same + //! IndexType. Is true by default if the BoxArray is empty. + bool ok () const; + + //! Return true if set of intersecting Boxes in BoxArray is null. + bool isDisjoint () const; + + //! Create a BoxList from this BoxArray. + BoxList boxList () const; + + //! True if the IntVect is within any of the Boxes in this BoxArray. + bool contains (const IntVect& v) const; + + //! \brief True if the Box is contained in this BoxArray(+ng). + //! The Box must also have the same IndexType as those in this BoxArray. + bool contains (const Box& b, bool assume_disjoint_ba = false, + const IntVect& ng = IntVect(0)) const; + + //! \brief True if all Boxes in ba are contained in this BoxArray(+ng). + bool contains (const BoxArray& ba, bool assume_disjoint_ba = false, + const IntVect& ng = IntVect(0)) const; + + //! Return smallest Box that contains all Boxes in this BoxArray. + Box minimalBox () const; + Box minimalBox (Long& npts_avg_box) const; + + //! \brief True if the Box intersects with this BoxArray(+ghostcells). + //! The Box must have the same IndexType as those in this BoxArray. + bool intersects (const Box& b, int ng = 0) const; + + bool intersects (const Box& b, const IntVect& ng) const; + + //! Return intersections of Box and BoxArray + std::vector< std::pair > intersections (const Box& bx) const; + + //! Return intersections of Box and BoxArray(+ghostcells). + std::vector< std::pair > intersections (const Box& bx, bool first_only, int ng) const; + + std::vector< std::pair > intersections (const Box& bx, bool first_only, const IntVect& ng) const; + + //! intersect Box and BoxArray, then store the result in isects + void intersections (const Box& bx, std::vector< std::pair >& isects) const; + + //! intersect Box and BoxArray(+ghostcells), then store the result in isects + void intersections (const Box& bx, std::vector< std::pair >& isects, + bool first_only, int ng) const; + + void intersections (const Box& bx, std::vector< std::pair >& isects, + bool first_only, const IntVect& ng) const; + + //! Return box - boxarray + BoxList complementIn (const Box& b) const; + void complementIn (BoxList& bl, const Box& b) const; + + //! Clear out the internal hash table used by intersections. + void clear_hash_bin () const; + + //! Change the BoxArray to one with no overlap and then simplify it (see the simplify function in BoxList). + void removeOverlap (bool simplify=true); + + //! whether two BoxArrays share the same data + static bool SameRefs (const BoxArray& lhs, const BoxArray& rhs) { return lhs.m_ref == rhs.m_ref; } + + struct RefID { + RefID () noexcept : data(nullptr) {} + explicit RefID (BARef* data_) noexcept : data(data_) {} + bool operator< (const RefID& rhs) const noexcept { return std::less()(data,rhs.data); } + bool operator== (const RefID& rhs) const noexcept { return data == rhs.data; } + bool operator!= (const RefID& rhs) const noexcept { return data != rhs.data; } + friend std::ostream& operator<< (std::ostream& os, const RefID& id); + private: + BARef* data; + }; + + //! Return a unique ID of the reference + RefID getRefID () const noexcept { return RefID { m_ref.get() }; } + + //! Return index type of this BoxArray + IndexType ixType () const noexcept { return m_bat.index_type(); } + + //! Return crse ratio of this BoxArray + IntVect crseRatio () const noexcept { return m_bat.coarsen_ratio(); } + + static void Initialize (); + static void Finalize (); + static bool initialized; + + //! Make ourselves unique. + void uniqify (); + + BoxList const& simplified_list () const; // For regular AMR grids only!!! + BoxArray simplified () const; // For regular AMR grids only!!! + +*/ + + ; +} diff --git a/src/Base/CMakeLists.txt b/src/Base/CMakeLists.txt index a7297c8f..c169d05b 100644 --- a/src/Base/CMakeLists.txt +++ b/src/Base/CMakeLists.txt @@ -2,7 +2,9 @@ target_sources(pyAMReX PRIVATE AMReX.cpp Box.cpp + BoxArray.cpp Dim3.cpp + DistributionMapping.cpp IntVect.cpp MultiFab.cpp ) diff --git a/src/Base/DistributionMapping.cpp b/src/Base/DistributionMapping.cpp new file mode 100644 index 00000000..5f57d0cc --- /dev/null +++ b/src/Base/DistributionMapping.cpp @@ -0,0 +1,79 @@ +/* Copyright 2021 The AMReX Community + * + * Authors: Axel Huebl + * License: BSD-3-Clause-LBNL + */ +#include +#include + +#include +#include +#include +#include + +#include + +namespace py = pybind11; +using namespace amrex; + + +void init_DistributionMapping(py::module &m) { + py::class_< DistributionMapping >(m, "DistributionMapping") + .def("__repr__", + [](DistributionMapping const & dm) { + std::stringstream s; + s << dm.size(); + return ""; + } + ) + + .def(py::init< >()) + .def(py::init< DistributionMapping const & >()) + //.def(py::init< DistributionMapping && >()) + //.def(py::init< DistributionMapping const &, DistributionMapping const & >()) + .def(py::init< Vector< int > const & >()) + //.def(py::init< Vector< int > && >()) + .def(py::init< BoxArray const & >(), + py::arg("boxes") + ) + .def(py::init< BoxArray const &, int >(), + py::arg("boxes"), py::arg("nprocs") + ) + + .def("define", + [](DistributionMapping & dm, BoxArray const & boxes) { + dm.define(boxes); + }, + py::arg("boxes") + ) + .def("define", + py::overload_cast< BoxArray const &, int >(&DistributionMapping::define), + py::arg("boxes"), py::arg("nprocs") + ) + .def("define", + py::overload_cast< Vector< int > const & >(&DistributionMapping::define)) + //.def("define", + // py::overload_cast< Vector< int > && >(&DistributionMapping::define)) + //! Length of the underlying processor map. + .def_property_readonly("size", &DistributionMapping::size) + .def_property_readonly("capacity", &DistributionMapping::capacity) + .def_property_readonly("empty", &DistributionMapping::empty) + + //! Number of references to this DistributionMapping + .def_property_readonly("link_count", &DistributionMapping::linkCount) + + /** + * \brief Returns a constant reference to the mapping of boxes in the + * underlying BoxArray to the CPU that holds the FAB on that Box. + * ProcessorMap()[i] is an integer in the interval [0, NCPU) where + * NCPU is the number of CPUs being used. + */ + .def("ProcessorMap", &DistributionMapping::ProcessorMap) + + //! Equivalent to ProcessorMap()[index]. + .def("__getitem__", + [](DistributionMapping const & dm, int index) -> int { + return dm[index]; + }) + ; +} diff --git a/src/pyAMReX.cpp b/src/pyAMReX.cpp index 32ee9ae3..2d5e0f2e 100644 --- a/src/pyAMReX.cpp +++ b/src/pyAMReX.cpp @@ -18,7 +18,9 @@ namespace py = pybind11; // forward declarations of exposed classes void init_AMReX(py::module&); void init_Box(py::module &); +void init_BoxArray(py::module &); void init_Dim3(py::module&); +void init_DistributionMapping(py::module&); void init_IntVect(py::module &); void init_MultiFab(py::module &); @@ -31,7 +33,10 @@ PYBIND11_MODULE(amrex_pybind, m) { .. autosummary:: :toctree: _generate Box + BoxArray + Dim3 IntVect + MultiFab )pbdoc"; // note: order from parent to child classes @@ -39,7 +44,9 @@ PYBIND11_MODULE(amrex_pybind, m) { init_Dim3(m); init_IntVect(m); init_Box(m); + init_BoxArray(m); init_MultiFab(m); + init_DistributionMapping(m); // API runtime version // note PEP-440 syntax: x.y.zaN but x.y.z.devN From 2ac5c81b71f7f6b6e5659190f1c34307368d30d9 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 7 Oct 2021 11:40:24 -0700 Subject: [PATCH 05/43] More Unittests Co-authored-by: Shreyas Ananthan --- tests/conftest.py | 15 ++++++++ tests/test_box.py | 69 ++++++++++++++++++++++++++++++++++ tests/test_multifab.py | 84 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 168 insertions(+) create mode 100644 tests/test_box.py create mode 100644 tests/test_multifab.py diff --git a/tests/conftest.py b/tests/conftest.py index b5ec37dc..275d1e32 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -11,3 +11,18 @@ def amrex_init(): amrex.initialize(["amrex.verbose=-1"]) yield amrex.finalize() + +@pytest.fixture(scope='module') +def boxarr(): + """BoxArray for MultiFab creation""" + #bx = amrex.Box.new((0, 0, 0), (63, 63, 63)) + bx = amrex.Box(amrex.IntVect(0, 0, 0), amrex.IntVect(63, 63, 63)) + ba = amrex.BoxArray(bx) + ba.max_size(32) + return ba + +@pytest.fixture(scope='module') +def distmap(boxarr): + """DistributionMapping for MultiFab creation""" + dm = amrex.DistributionMapping(boxarr) + return dm diff --git a/tests/test_box.py b/tests/test_box.py new file mode 100644 index 00000000..0eb373fa --- /dev/null +++ b/tests/test_box.py @@ -0,0 +1,69 @@ +# -*- coding: utf-8 -*- + +import pytest +import numpy as np +import amrex + +@pytest.fixture +def box(): + #return amrex.Box((0, 0, 0), (127, 127, 127)) + return amrex.Box(amrex.IntVect(0, 0, 0), amrex.IntVect(127, 127, 127)) + +#def test_num_pts(box): +# np.testing.assert_allclose(box.lo_vect, [0, 0, 0]) +# np.testing.assert_allclose(box.hi_vect, [127, 127, 127]) +# assert(box.num_pts == 2**21) +# assert(box.volume == 2**21) + +#def test_grow(box): +# """box.grow""" +# bx = box.grow(3) +# np.testing.assert_allclose(bx.lo_vect, [-3, -3, -3]) +# np.testing.assert_allclose(bx.hi_vect, [130, 130, 130]) +# assert(bx.num_pts == (134**3)) + +#def test_convert(box): +# """Conversion to node""" +# bx = box.convert(amrex.CellIndex.NODE, amrex.CellIndex.NODE, amrex.CellIndex.NODE) +# assert(bx.num_pts == 129**3) +# assert(bx.volume == 128**3) + +# bx = box.convert(amrex.CellIndex.NODE, amrex.CellIndex.CELL, amrex.CellIndex.CELL) +# np.testing.assert_allclose(bx.hi_vect, [128, 127, 127]) +# assert(bx.num_pts == 129 * 128 * 128) +# assert(bx.volume == 128**3) + +@pytest.mark.parametrize("dir", [-1, 0, 1, 2]) +def test_surrounding_nodes(box, dir): + """Surrounding nodes""" + nx = np.array(box.hi_vect) + bx = box.surrounding_nodes(dir=dir) + + if dir < 0: + assert(bx.num_pts == 129**3) + assert(bx.volume == 128**3) + nx += 1 + np.testing.assert_allclose(bx.hi_vect, nx) + else: + assert(bx.num_pts == 129 * 128 * 128) + assert(bx.volume == 128**3) + nx[dir] += 1 + np.testing.assert_allclose(bx.hi_vect, nx) + +@pytest.mark.parametrize("dir", [-1, 0, 1, 2]) +def test_enclosed_cells(box, dir): + """Enclosed cells""" + bxn = box.convert(amrex.CellIndex.NODE, amrex.CellIndex.NODE, amrex.CellIndex.NODE) + nx = np.array(bxn.hi_vect) + bx = bxn.enclosed_cells(dir) + + if dir < 0: + assert(bx.num_pts == 128**3) + assert(bx.volume == 128**3) + nx -= 1 + np.testing.assert_allclose(bx.hi_vect, nx) + else: + assert(bx.num_pts == 129 * 129 * 128) + assert(bx.volume == 128**3) + nx[dir] -= 1 + np.testing.assert_allclose(bx.hi_vect, nx) diff --git a/tests/test_multifab.py b/tests/test_multifab.py new file mode 100644 index 00000000..e80ba8f2 --- /dev/null +++ b/tests/test_multifab.py @@ -0,0 +1,84 @@ +# -*- coding: utf-8 -*- + +import itertools +import numpy as np +import pytest +import amrex + + +#@pytest.fixture(params=list(itertools.product([1,3],[0,1]))) +def mfab1(boxarr, distmap, request): + """MultiFab for tests""" + num_components = 3 #request.param[0] + num_ghost = 1 #request.param[1] + mfab = amrex.MultiFab(boxarr, distmap, num_components, num_ghost) + mfab.set_val(0.0, 0, num_components) + return mfab # issue here when used as fixture + +@pytest.mark.parametrize("nghost", [0, 1]) +def test_mfab_loop(boxarr, distmap, nghost): + mfab = amrex.MultiFab(boxarr, distmap, 3, nghost) + # Looping + for mfi in mfab: + bx = mfi.tilebox() + marr = mfab.array(mfi) # Array4: add __array_interface__ here + + for i, j, k in bx: + mar[i, j, k, 0] = 10.0 * i + mar[i, j, k, 1] = 10.0 * j + mar[i, j, k, 2] = 10.0 * k + print(i,j,k) + +def test_mfab_simple(boxarr, distmap, request): + mfab = mfab1(boxarr, distmap, request) + assert(mfab.is_cell_centered) + assert(all(not mfab.is_nodal(i) for i in [-1, 0, 1, 2])) + + for i in range(mfab.num_comp): + mfab.set_val(-10 * (i + 1), i, 1) + mfab.abs(0, mfab.num_comp) + for i in range(mfab.num_comp): + assert(mfab.max(i) == (10 * (i + 1))) # Assert: None == 10 for i=0 + assert(mfab.min(i) == (10 * (i + 1))) + + mfab.plus(20.0, 0, mfab.num_comp) + for i in range(mfab.num_comp): + np.testing.assert_allclose(mfab.max(i), 20.0 + (10 * (i + 1))) + np.testing.assert_allclose(mfab.min(i), 20.0 + (10 * (i + 1))) + + mfab.mult(10.0, 0, mfab.num_comp) + for i in range(mfab.num_comp): + np.testing.assert_allclose(mfab.max(i), 10.0 * (20.0 + (10 * (i + 1)))) + np.testing.assert_allclose(mfab.min(i), 10.0 * (20.0 + (10 * (i + 1)))) + mfab.invert(10.0, 0, mfab.num_comp) + for i in range(mfab.num_comp): + np.testing.assert_allclose(mfab.max(i), 1.0 / (20.0 + (10 * (i + 1)))) + np.testing.assert_allclose(mfab.min(i), 1.0 / (20.0 + (10 * (i + 1)))) +''' +@pytest.mark.parametrize("nghost", [0, 1]) +def test_mfab_ops(boxarr, distmap, nghost): + src = amrex.MultiFab(boxarr, distmap, 3, nghost) + dst = amrex.MultiFab(boxarr, distmap, 1, nghost) + + src.set_val(10.0, 0, 1) + src.set_val(20.0, 1, 1) + src.set_val(30.0, 2, 1) + dst.set_val(0.0, 0, 1) + + dst.add(src, 2, 0, 1, nghost) + dst.subtract(src, 1, 0, 1, nghost) + dst.multiply(src, 0, 0, 1, nghost) + dst.divide(src, 1, 0, 1, nghost) + np.testing.assert_allclose(dst.min(0), 5.0) + np.testing.assert_allclose(dst.max(0), 5.0) + + dst.xpay(2.0, src, 0, 0, 1, nghost) + dst.saxpy(2.0, src, 1, 0, 1, nghost) + np.testing.assert_allclose(dst.min(0), 60.0) + np.testing.assert_allclose(dst.max(0), 60.0) + + dst.lin_comb(6.0, src, 1, + 1.0, src, 2, 0, 1, nghost) + np.testing.assert_allclose(dst.min(0), 150.0) + np.testing.assert_allclose(dst.max(0), 150.0) +''' From 020caf95800c2e93f21872d36e3173b48c57e43b Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 7 Oct 2021 18:57:23 -0700 Subject: [PATCH 06/43] MultiFab: Fix tests Forgot to return in custom lambda methods --- src/Base/MultiFab.cpp | 10 +++++++-- tests/conftest.py | 10 +++++++++ tests/test_multifab.py | 51 ++++++++++++++++++++++-------------------- 3 files changed, 45 insertions(+), 26 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index b6ce8a03..9ca493ba 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -111,14 +111,14 @@ void init_MultiFab(py::module &m) { /* sizes, etc. */ .def("min", - [](MultiFab const & mf, int comp) { mf.min(comp); }) + [](MultiFab const & mf, int comp) { return mf.min(comp); }) .def("min", py::overload_cast< int, int, bool >(&MultiFab::min, py::const_)) .def("min", py::overload_cast< Box const &, int, int, bool >(&MultiFab::min, py::const_)) .def("max", - [](MultiFab const & mf, int comp) { mf.max(comp); }) + [](MultiFab const & mf, int comp) { return mf.max(comp); }) .def("max", py::overload_cast< int, int, bool >(&MultiFab::max, py::const_)) .def("max", @@ -150,6 +150,8 @@ void init_MultiFab(py::module &m) { .def("abs", py::overload_cast< int, int, int >(&MultiFab::abs)) .def("plus", py::overload_cast< Real, int >(&MultiFab::plus)) + .def("plus", [](MultiFab & mf, Real val, int comp, int num_comp) { + mf.plus(val, comp, num_comp); }) .def("plus", py::overload_cast< Real, int, int, int >(&MultiFab::plus)) .def("plus", py::overload_cast< Real, Box const &, int >(&MultiFab::plus)) .def("plus", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::plus)) @@ -161,11 +163,15 @@ void init_MultiFab(py::module &m) { .def("divi", py::overload_cast< MultiFab const &, int, int, int >(&MultiFab::divide)) .def("mult", py::overload_cast< Real, int >(&MultiFab::mult)) + .def("mult", [](MultiFab & mf, Real val, int comp, int num_comp) { + mf.mult(val, comp, num_comp); }) .def("mult", py::overload_cast< Real, int, int, int >(&MultiFab::mult)) .def("mult", py::overload_cast< Real, Box const &, int >(&MultiFab::mult)) .def("mult", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::mult)) .def("invert", py::overload_cast< Real, int >(&MultiFab::invert)) + .def("invert", [](MultiFab & mf, Real val, int comp, int num_comp) { + mf.invert(val, comp, num_comp); }) .def("invert", py::overload_cast< Real, int, int, int >(&MultiFab::invert)) .def("invert", py::overload_cast< Real, Box const &, int >(&MultiFab::invert)) .def("invert", py::overload_cast< Real, Box const &, int, int, int >(&MultiFab::invert)) diff --git a/tests/conftest.py b/tests/conftest.py index 275d1e32..5f48aa68 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,5 +1,6 @@ # -*- coding: utf-8 -*- +import itertools import pytest import amrex @@ -26,3 +27,12 @@ def distmap(boxarr): """DistributionMapping for MultiFab creation""" dm = amrex.DistributionMapping(boxarr) return dm + +@pytest.fixture(params=list(itertools.product([1,3],[0,1]))) +def mfab(boxarr, distmap, request): + """MultiFab for tests""" + num_components = request.param[0] + num_ghost = request.param[1] + mfab = amrex.MultiFab(boxarr, distmap, num_components, num_ghost) + mfab.set_val(0.0, 0, num_components) + return mfab diff --git a/tests/test_multifab.py b/tests/test_multifab.py index e80ba8f2..18c1838a 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -1,25 +1,14 @@ # -*- coding: utf-8 -*- -import itertools import numpy as np import pytest import amrex -#@pytest.fixture(params=list(itertools.product([1,3],[0,1]))) -def mfab1(boxarr, distmap, request): - """MultiFab for tests""" - num_components = 3 #request.param[0] - num_ghost = 1 #request.param[1] - mfab = amrex.MultiFab(boxarr, distmap, num_components, num_ghost) - mfab.set_val(0.0, 0, num_components) - return mfab # issue here when used as fixture - @pytest.mark.parametrize("nghost", [0, 1]) -def test_mfab_loop(boxarr, distmap, nghost): - mfab = amrex.MultiFab(boxarr, distmap, 3, nghost) - # Looping +def test_mfab_loop(mfab, nghost): for mfi in mfab: + print(mfi) bx = mfi.tilebox() marr = mfab.array(mfi) # Array4: add __array_interface__ here @@ -28,9 +17,10 @@ def test_mfab_loop(boxarr, distmap, nghost): mar[i, j, k, 1] = 10.0 * j mar[i, j, k, 2] = 10.0 * k print(i,j,k) + assert(mar[1, 2, 3, 1] == 42) + -def test_mfab_simple(boxarr, distmap, request): - mfab = mfab1(boxarr, distmap, request) +def test_mfab_simple(mfab): assert(mfab.is_cell_centered) assert(all(not mfab.is_nodal(i) for i in [-1, 0, 1, 2])) @@ -54,7 +44,8 @@ def test_mfab_simple(boxarr, distmap, request): for i in range(mfab.num_comp): np.testing.assert_allclose(mfab.max(i), 1.0 / (20.0 + (10 * (i + 1)))) np.testing.assert_allclose(mfab.min(i), 1.0 / (20.0 + (10 * (i + 1)))) -''' + + @pytest.mark.parametrize("nghost", [0, 1]) def test_mfab_ops(boxarr, distmap, nghost): src = amrex.MultiFab(boxarr, distmap, 3, nghost) @@ -65,20 +56,32 @@ def test_mfab_ops(boxarr, distmap, nghost): src.set_val(30.0, 2, 1) dst.set_val(0.0, 0, 1) - dst.add(src, 2, 0, 1, nghost) - dst.subtract(src, 1, 0, 1, nghost) - dst.multiply(src, 0, 0, 1, nghost) - dst.divide(src, 1, 0, 1, nghost) + #dst.add(src, 2, 0, 1, nghost) + #dst.subtract(src, 1, 0, 1, nghost) + #dst.multiply(src, 0, 0, 1, nghost) + #dst.divide(src, 1, 0, 1, nghost) + + dst.add(dst, src, 2, 0, 1, nghost) + dst.subtract(dst, src, 1, 0, 1, nghost) + dst.multiply(dst, src, 0, 0, 1, nghost) + dst.divide(dst, src, 1, 0, 1, nghost) + + print(dst.min(0)) np.testing.assert_allclose(dst.min(0), 5.0) np.testing.assert_allclose(dst.max(0), 5.0) - dst.xpay(2.0, src, 0, 0, 1, nghost) - dst.saxpy(2.0, src, 1, 0, 1, nghost) + #dst.xpay(2.0, src, 0, 0, 1, nghost) + #dst.saxpy(2.0, src, 1, 0, 1, nghost) + dst.xpay(dst, 2.0, src, 0, 0, 1, nghost) + dst.saxpy(dst, 2.0, src, 1, 0, 1, nghost) np.testing.assert_allclose(dst.min(0), 60.0) np.testing.assert_allclose(dst.max(0), 60.0) - dst.lin_comb(6.0, src, 1, + #dst.lin_comb(6.0, src, 1, + # 1.0, src, 2, 0, 1, nghost) + dst.lin_comb(dst, + 6.0, src, 1, 1.0, src, 2, 0, 1, nghost) np.testing.assert_allclose(dst.min(0), 150.0) np.testing.assert_allclose(dst.max(0), 150.0) -''' + From e1f78099bcf289cab57b83dd78dbe255276c193f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 7 Oct 2021 19:55:57 -0700 Subject: [PATCH 07/43] MFIter: Close Still fails with double free, probably because MFIter survives MultiFab or so --- src/Base/MultiFab.cpp | 47 ++++++++++++++++++++++++++++++++++++++++-- tests/test_multifab.py | 17 +++++++-------- 2 files changed, 53 insertions(+), 11 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 9ca493ba..c2eca89a 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -41,7 +41,7 @@ namespace { class ValidSentinel {}; bool operator==(MFIterWrapper const& it, ValidSentinel const&) { - return (*it)->isValid(); + return !(*it)->isValid(); } } @@ -50,6 +50,49 @@ void init_MultiFab(py::module &m) { py::class_< FArrayBox >(m, "FArrayBox"); py::class_< FabArray, FabArrayBase >(m, "FabArray_FArrayBox"); + //py::class_< MFIterWrapper >(m, "MFIterWrapper"); + py::class_< MFIter >(m, "MFIter") + .def("__repr__", + [](MFIter const & mfi) { + std::string r = ""); + return r; + } + ) + .def(py::init< FabArrayBase const & >()) + .def(py::init< FabArrayBase const &, MFItInfo const & >()) + + .def(py::init< MultiFab const & >()) + .def(py::init< MultiFab const &, MFItInfo const & >()) + + //.def(py::init< iMultiFab const & >()) + //.def(py::init< iMultiFab const &, MFItInfo const & >()) + + .def("tilebox", py::overload_cast< >(&MFIter::tilebox, py::const_)) + .def("tilebox", py::overload_cast< IntVect const & >(&MFIter::tilebox, py::const_)) + .def("tilebox", py::overload_cast< IntVect const &, IntVect const & >(&MFIter::tilebox, py::const_)) + + /* + Box nodaltilebox() + Box nodaltilebox(int dir) + Box growntilebox() + Box growntilebox(const IntVect&) + Box grownnodaltilebox() + Box grownnodaltilebox(int dir) + Box grownnodaltilebox(int dir, int ng) + Box grownnodaltilebox(int dir, const IntVect&) + + Box validbox() + Box fabbox() + + void operator++() + bint isValid() + int index() + int length() + */ + ; + py::class_< MultiFab /*, FabArray*/ >(m, "MultiFab") .def("__repr__", [](MultiFab const & mf) { @@ -246,7 +289,7 @@ void init_MultiFab(py::module &m) { /* data access in Box index space */ .def("__iter__", - [](const MultiFab& mf) { + [](MultiFab& mf) { return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); }, // Essential: keep object alive while iterator exists diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 18c1838a..1cc423ab 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -9,15 +9,14 @@ def test_mfab_loop(mfab, nghost): for mfi in mfab: print(mfi) - bx = mfi.tilebox() - marr = mfab.array(mfi) # Array4: add __array_interface__ here - - for i, j, k in bx: - mar[i, j, k, 0] = 10.0 * i - mar[i, j, k, 1] = 10.0 * j - mar[i, j, k, 2] = 10.0 * k - print(i,j,k) - assert(mar[1, 2, 3, 1] == 42) + #bx = mfi.tilebox() + #marr = mfab.array(mfi) # Array4: add __array_interface__ here + + #for i, j, k in bx: + # mar[i, j, k, 0] = 10.0 * i + # mar[i, j, k, 1] = 10.0 * j + # mar[i, j, k, 2] = 10.0 * k + # print(i,j,k) def test_mfab_simple(mfab): From f9004024f721ecd4eb1f52f3f2438e63ecd65d87 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 8 Oct 2021 00:26:32 -0700 Subject: [PATCH 08/43] Implement: MFIter --- src/Base/MultiFab.cpp | 44 +++++++++++++++++------------------------- tests/test_multifab.py | 2 +- 2 files changed, 19 insertions(+), 27 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index c2eca89a..3f8067e7 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -19,31 +19,6 @@ namespace py = pybind11; using namespace amrex; -namespace { - /** STL-like iteration over amrex::MFIter - * - * The amrex::MFIter interface is currently a bit too tricky to implement - * std::begin() and std::end() safely with OpenMP threading. - */ - class MFIterWrapper { - std::shared_ptr mfi; - public: - explicit MFIterWrapper(const MultiFab& mf) { - // For tiling support (OpenMP/thread pools) later on: - // MFIter mfi(mf, TilingIfNotGPU()); - mfi = std::make_shared(mf); - } - std::shared_ptr operator*() { return mfi; } - std::shared_ptr operator*() const { return mfi; } - MFIterWrapper& operator++() { ++(*mfi); return *this; } - }; - - class ValidSentinel {}; - - bool operator==(MFIterWrapper const& it, ValidSentinel const&) { - return !(*it)->isValid(); - } -} void init_MultiFab(py::module &m) { py::class_< FabArrayBase >(m, "FabArrayBase"); @@ -69,6 +44,22 @@ void init_MultiFab(py::module &m) { //.def(py::init< iMultiFab const & >()) //.def(py::init< iMultiFab const &, MFItInfo const & >()) + //.def("__iter__", + // [](MFIter & mfi) -> MFIter & { + // return mfi; + // }, + // py::return_value_policy::reference + //) + .def("__next__", + [](MFIter & mfi) -> MFIter & { + ++mfi; + if( !mfi.isValid() ) + throw py::stop_iteration(); + return mfi; + }, + py::return_value_policy::reference + ) + .def("tilebox", py::overload_cast< >(&MFIter::tilebox, py::const_)) .def("tilebox", py::overload_cast< IntVect const & >(&MFIter::tilebox, py::const_)) .def("tilebox", py::overload_cast< IntVect const &, IntVect const & >(&MFIter::tilebox, py::const_)) @@ -290,7 +281,8 @@ void init_MultiFab(py::module &m) { /* data access in Box index space */ .def("__iter__", [](MultiFab& mf) { - return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); + //return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); + return MFIter(mf); }, // Essential: keep object alive while iterator exists py::keep_alive<0, 1>() diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 1cc423ab..1ab636c0 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -9,7 +9,7 @@ def test_mfab_loop(mfab, nghost): for mfi in mfab: print(mfi) - #bx = mfi.tilebox() + bx = mfi.tilebox() #marr = mfab.array(mfi) # Array4: add __array_interface__ here #for i, j, k in bx: From 77e1cf6024208ed3b0048cfeaefdcdcac6dc657c Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 8 Oct 2021 00:52:22 -0700 Subject: [PATCH 09/43] Box Updates --- src/Base/Box.cpp | 23 +++++++++++++++++++---- tests/test_box.py | 2 ++ 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/Base/Box.cpp b/src/Base/Box.cpp index 5926e0f2..022a3e62 100644 --- a/src/Base/Box.cpp +++ b/src/Base/Box.cpp @@ -17,6 +17,8 @@ using namespace amrex; void init_Box(py::module &m) { + py::class_< Direction >(m, "Direction"); + py::class_< Box >(m, "Box") .def("__repr__", [](Box const & b) { @@ -57,8 +59,8 @@ void init_Box(py::module &m) { // loVect3d // hiVect3d - // loVect - // hiVect + .def_property_readonly("lo_vect", &Box::loVect) + .def_property_readonly("hi_vect", &Box::hiVect) .def("contains", py::overload_cast< IntVect const & >(&Box::contains, py::const_)) @@ -76,8 +78,21 @@ void init_Box(py::module &m) { // setRange // shift // shiftHalf - // convert - // surroundingNodes + + .def("convert", + py::overload_cast< IndexType >(&Box::convert)) + .def("convert", + py::overload_cast< IntVect const & >(&Box::convert)) + + //.def("surrounding_nodes", + // py::overload_cast< >(&Box::surroundingNodes)) + //.def("surrounding_nodes", + // py::overload_cast< int >(&Box::surroundingNodes), + // py::arg("dir")) + //.def("surrounding_nodes", + // py::overload_cast< Direction >(&Box::surroundingNodes), + // py::arg("d")) + // enclosedCells // minBox // chop diff --git a/tests/test_box.py b/tests/test_box.py index 0eb373fa..05b854c6 100644 --- a/tests/test_box.py +++ b/tests/test_box.py @@ -33,6 +33,7 @@ def box(): # assert(bx.num_pts == 129 * 128 * 128) # assert(bx.volume == 128**3) +''' @pytest.mark.parametrize("dir", [-1, 0, 1, 2]) def test_surrounding_nodes(box, dir): """Surrounding nodes""" @@ -67,3 +68,4 @@ def test_enclosed_cells(box, dir): assert(bx.volume == 128**3) nx[dir] -= 1 np.testing.assert_allclose(bx.hi_vect, nx) +''' From 766b30022562d216fd01cdcabdcd9b36673cc6f9 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 18:32:39 -0800 Subject: [PATCH 10/43] MultiFab: Compile Issues Try again with a newer pybind11 version later? --- src/Base/MultiFab.cpp | 45 +++++++++++++++++++++++++++++------------- tests/test_multifab.py | 2 +- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 3f8067e7..3aefbd01 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -21,8 +21,18 @@ using namespace amrex; void init_MultiFab(py::module &m) { - py::class_< FabArrayBase >(m, "FabArrayBase"); + py::class_< FabArrayBase >(m, "FabArrayBase") + .def_property_readonly("is_all_cell_centered", &FabArrayBase::is_cell_centered) + .def_property_readonly("is_all_nodal", + py::overload_cast< >(&FabArrayBase::is_nodal, py::const_)) + .def("is_nodal", + py::overload_cast< int >(&FabArrayBase::is_nodal, py::const_)) + + .def_property_readonly("num_comp", &FabArrayBase::nComp) + ; + py::class_< FArrayBox >(m, "FArrayBox"); + py::class_< FabArray, FabArrayBase >(m, "FabArray_FArrayBox"); //py::class_< MFIterWrapper >(m, "MFIterWrapper"); @@ -84,7 +94,7 @@ void init_MultiFab(py::module &m) { */ ; - py::class_< MultiFab /*, FabArray*/ >(m, "MultiFab") + py::class_< MultiFab, FabArray >(m, "MultiFab") .def("__repr__", [](MultiFab const & mf) { return "(&MultiFab::define)) /* setters */ + //.def("set_val", + // py::overload_cast< Real >(&MultiFab::setVal)) .def("set_val", - py::overload_cast< Real >(&MultiFab::setVal)) + [](MultiFab & mf, Real val) { mf.setVal(val); } + ) .def("set_val", [](MultiFab & mf, Real val, int comp, int num_comp) { mf.setVal(val, comp, num_comp); } ) + //.def("set_val", + // py::overload_cast< Real, int, int, int >(&MultiFab::setVal)) .def("set_val", - py::overload_cast< Real, int, int, int >(&MultiFab::setVal)) + [](MultiFab & mf, Real val, int comp, int ncomp, int nghost) { + mf.setVal(val, comp, ncomp, nghost); + } + ) + //.def("set_val", + // py::overload_cast< Real, int, int, IntVect const & >(&MultiFab::setVal)) .def("set_val", - py::overload_cast< Real, int, int, IntVect const & >(&MultiFab::setVal)) - - .def("is_cell_centered", &MultiFab::is_cell_centered) - .def("is_nodal", - py::overload_cast< >(&MultiFab::is_nodal, py::const_)) - .def("is_nodal", - py::overload_cast< int >(&MultiFab::is_nodal, py::const_)) - - .def_property_readonly("num_comp", &MultiFab::nComp) + [](MultiFab & mf, Real val, int comp, int ncomp, IntVect const & nghost) { + mf.setVal(val, comp, ncomp, nghost); + } + ) /* sizes, etc. */ .def("min", @@ -181,7 +196,9 @@ void init_MultiFab(py::module &m) { .def("abs", [](MultiFab & mf, int comp, int num_comp) { mf.abs(comp, num_comp); }) - .def("abs", py::overload_cast< int, int, int >(&MultiFab::abs)) + //.def("abs", py::overload_cast< int, int, int >(&MultiFab::abs)) + .def("abs", + [](MultiFab & mf, int comp, int num_comp, int nghost) { mf.abs(comp, num_comp, nghost); }) .def("plus", py::overload_cast< Real, int >(&MultiFab::plus)) .def("plus", [](MultiFab & mf, Real val, int comp, int num_comp) { diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 1ab636c0..f0167788 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -20,7 +20,7 @@ def test_mfab_loop(mfab, nghost): def test_mfab_simple(mfab): - assert(mfab.is_cell_centered) + assert(mfab.is_all_cell_centered) assert(all(not mfab.is_nodal(i) for i in [-1, 0, 1, 2])) for i in range(mfab.num_comp): From 632cff8e27c3cd69d4d06b6898a5fab366f146c2 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 15 Feb 2021 01:27:33 -0800 Subject: [PATCH 11/43] [Draft] Array4: __array_interface__ Implement the `__array_interface__` in Array for. --- src/Base/Array4.cpp | 160 ++++++++++++++++++++++++++++++++++++++++ src/Base/CMakeLists.txt | 1 + src/pyAMReX.cpp | 2 + tests/test_array4.py | 86 +++++++++++++++++++++ 4 files changed, 249 insertions(+) create mode 100644 src/Base/Array4.cpp create mode 100644 tests/test_array4.py diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp new file mode 100644 index 00000000..caf53369 --- /dev/null +++ b/src/Base/Array4.cpp @@ -0,0 +1,160 @@ +/* Copyright 2021 The AMReX Community + * + * Authors: Axel Huebl + * License: BSD-3-Clause-LBNL + */ +#include +#include +#include + +#include +#include +#include + +#include +#include + +namespace py = pybind11; +using namespace amrex; + + +template< typename T > +void make_Array4(py::module &m, std::string typestr) +{ + // dispatch simpler via: py::format_descriptor::format() naming + auto const array_name = std::string("Array4_").append(typestr); + py::class_< Array4 >(m, array_name.c_str(), py::buffer_protocol()) + .def("__repr__", + [](Array4 const & a4) { + std::stringstream s; + s << a4.size(); + return ""; + } + ) +#if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK) + .def("index_assert", &Array4::index_assert) +#endif + + .def_property_readonly("size", &Array4::size) + .def_property_readonly("nComp", &Array4::nComp) + + .def(py::init< >()) + .def(py::init< Array4 const & >()) + .def(py::init< Array4 const &, int >()) + .def(py::init< Array4 const &, int, int >()) + //.def(py::init< T*, Dim3 const &, Dim3 const &, int >()) + + .def(py::init([](py::array_t & arr) { + py::buffer_info buf = arr.request(); + + auto a4 = std::make_unique< Array4 >(); + a4.get()->p = (T*)buf.ptr; + a4.get()->begin = Dim3{0, 0, 0}; + // TODO: likely C->F index conversion here + // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; + a4.get()->end.x = (int)buf.shape.at(0); + a4.get()->end.y = (int)buf.shape.at(1); + a4.get()->end.z = (int)buf.shape.at(2); + a4.get()->ncomp = 1; + // buffer protocol strides are in bytes, AMReX strides are elements + a4.get()->jstride = (int)buf.strides.at(0) / sizeof(T); + a4.get()->kstride = (int)buf.strides.at(1) / sizeof(T); + a4.get()->nstride = (int)buf.strides.at(2) * (int)buf.shape.at(2) / sizeof(T); + return a4; + })) + + + .def_property_readonly("__array_interface__", [](Array4 const & a4) { + auto d = py::dict(); + auto const len = length(a4); + // TODO: likely F->C index conversion here + // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; + auto shape = py::make_tuple( // Buffer dimensions + len.x < 0 ? 0 : len.x, + len.y < 0 ? 0 : len.y, + len.z < 0 ? 0 : len.z//, // zero-size shall not have negative dimension + //a4.ncomp + ); + // buffer protocol strides are in bytes, AMReX strides are elements + auto const strides = py::make_tuple( // Strides (in bytes) for each index + sizeof(T) * a4.jstride, + sizeof(T) * a4.kstride, + sizeof(T)//, + //sizeof(T) * a4.nstride + ); + d["data"] = py::make_tuple(long(a4.dataPtr()), false); + d["typestr"] = py::format_descriptor::format(); + d["shape"] = shape; + d["strides"] = strides; + // d["strides"] = py::none(); + d["version"] = 3; + return d; + }) + + // not sure if useful to have this implemented on top +/* + .def_buffer([](Array4 & a4) -> py::buffer_info { + auto const len = length(a4); + // TODO: likely F->C index conversion here + // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; + auto shape = { // Buffer dimensions + len.x < 0 ? 0 : len.x, + len.y < 0 ? 0 : len.y, + len.z < 0 ? 0 : len.z//, // zero-size shall not have negative dimension + //a4.ncomp + }; + // buffer protocol strides are in bytes, AMReX strides are elements + auto const strides = { // Strides (in bytes) for each index + sizeof(T) * a4.jstride, + sizeof(T) * a4.kstride, + sizeof(T)//, + //sizeof(T) * a4.nstride + }; + return py::buffer_info( + a4.dataPtr(), + shape, + strides + ); + }) +*/ + ; +} + +void init_Array4(py::module &m) { + make_Array4< float >(m, "float"); + make_Array4< double >(m, "double"); + make_Array4< long double >(m, "longdouble"); + + make_Array4< short >(m, "short"); + make_Array4< int >(m, "int"); + make_Array4< long >(m, "long"); + make_Array4< long long >(m, "longlong"); + + make_Array4< unsigned short >(m, "ushort"); + make_Array4< unsigned int >(m, "uint"); + make_Array4< unsigned long >(m, "ulong"); + make_Array4< unsigned long long >(m, "ulonglong"); + + // std::complex< float|double|long double> ? + + /* + py::class_< PolymorphicArray4, Array4 >(m, "PolymorphicArray4") + .def("__repr__", + [](PolymorphicArray4 const & pa4) { + std::stringstream s; + s << pa4.size(); + return ""; + } + ) + ; + */ + + // free standing C++ functions: + /* + contains + lbound + ubound + length + makePolymorphic + */ +} diff --git a/src/Base/CMakeLists.txt b/src/Base/CMakeLists.txt index c169d05b..2a82c2d3 100644 --- a/src/Base/CMakeLists.txt +++ b/src/Base/CMakeLists.txt @@ -1,6 +1,7 @@ target_sources(pyAMReX PRIVATE AMReX.cpp + Array4.cpp Box.cpp BoxArray.cpp Dim3.cpp diff --git a/src/pyAMReX.cpp b/src/pyAMReX.cpp index 2d5e0f2e..9a7987c5 100644 --- a/src/pyAMReX.cpp +++ b/src/pyAMReX.cpp @@ -17,6 +17,7 @@ namespace py = pybind11; // forward declarations of exposed classes void init_AMReX(py::module&); +void init_Array4(py::module&); void init_Box(py::module &); void init_BoxArray(py::module &); void init_Dim3(py::module&); @@ -43,6 +44,7 @@ PYBIND11_MODULE(amrex_pybind, m) { init_AMReX(m); init_Dim3(m); init_IntVect(m); + init_Array4(m); init_Box(m); init_BoxArray(m); init_MultiFab(m); diff --git a/tests/test_array4.py b/tests/test_array4.py new file mode 100644 index 00000000..0c8f0a42 --- /dev/null +++ b/tests/test_array4.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- + +import pytest +import numpy as np +import amrex + +def test_array4_empty(): + empty = amrex.Array4_double() + + # Check properties + assert(empty.size == 0) + assert(empty.nComp == 0) + + # assign empty + emptyc = amrex.Array4_double(empty) + # Check properties + assert(emptyc.size == 0) + assert(emptyc.nComp == 0) + +def test_array4(): + # from numpy (also a non-owning view) + x = np.ones((2, 3, 4,)) + arr = amrex.Array4_double(x) + + x[1, 1, 1] = 42 + # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable + # assert(arr[1, 1, 1] == 42) + + # copy to numpy + c_arr2np = np.array(arr, copy=True) + assert(c_arr2np.ndim == 3) + #assert(c_arr2np.ndim == 4) + assert(c_arr2np.dtype == np.dtype("double")) + np.testing.assert_array_equal(x, c_arr2np) + assert(c_arr2np[1, 1, 1] == 42) + + # view to numpy + v_arr2np = np.array(arr, copy=False) + assert(v_arr2np.ndim == 3) + #assert(c_arr2np.ndim == 4) + assert(v_arr2np.dtype == np.dtype("double")) + np.testing.assert_array_equal(x, v_arr2np) + assert(v_arr2np[1, 1, 1] == 42) + + # change original buffer once more + x[1, 1, 1] = 43 + assert(v_arr2np[1, 1, 1] == 43) + + # copy array4 (view) + c_arr = amrex.Array4_double(arr) + v_carr2np = np.array(c_arr, copy=False) + x[1, 1, 1] = 44 + assert(v_carr2np[1, 1, 1] == 44) + + # from cupy + + # to numpy + + # to cupy + + return + + # Check indexing + assert(obj[0] == 1) + assert(obj[1] == 2) + assert(obj[2] == 3) + assert(obj[-1] == 3) + assert(obj[-2] == 2) + assert(obj[-3] == 1) + with pytest.raises(IndexError): + obj[-4] + with pytest.raises(IndexError): + obj[3] + + # Check assignment + obj[0] = 2 + obj[1] = 3 + obj[2] = 4 + assert(obj[0] == 2) + assert(obj[1] == 3) + assert(obj[2] == 4) + +#def test_iv_conversions(): +# obj = amrex.IntVect.max_vector().numpy() +# assert(isinstance(obj, np.ndarray)) +# assert(obj.dtype == np.int32) From 0c32d1025747f61375d55c1587a10b4457e336b2 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 21:19:55 -0800 Subject: [PATCH 12/43] Array4 Updates (incl. ncomp) --- src/Base/Array4.cpp | 41 +++++++++++++++++++++++++++-------------- tests/test_array4.py | 11 +++++------ 2 files changed, 32 insertions(+), 20 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index caf53369..da8e8f72 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -37,6 +37,7 @@ void make_Array4(py::module &m, std::string typestr) .def_property_readonly("size", &Array4::size) .def_property_readonly("nComp", &Array4::nComp) + .def_property_readonly("num_comp", &Array4::nComp) .def(py::init< >()) .def(py::init< Array4 const & >()) @@ -44,6 +45,8 @@ void make_Array4(py::module &m, std::string typestr) .def(py::init< Array4 const &, int, int >()) //.def(py::init< T*, Dim3 const &, Dim3 const &, int >()) + /* init from a numpy or other buffer protocol array: non-owning view + */ .def(py::init([](py::array_t & arr) { py::buffer_info buf = arr.request(); @@ -63,24 +66,25 @@ void make_Array4(py::module &m, std::string typestr) return a4; })) - .def_property_readonly("__array_interface__", [](Array4 const & a4) { auto d = py::dict(); auto const len = length(a4); // TODO: likely F->C index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; + py::print("ncomp"); + py::print(a4.ncomp); auto shape = py::make_tuple( // Buffer dimensions len.x < 0 ? 0 : len.x, len.y < 0 ? 0 : len.y, - len.z < 0 ? 0 : len.z//, // zero-size shall not have negative dimension - //a4.ncomp + len.z < 0 ? 0 : len.z, // zero-size shall not have negative dimension + a4.ncomp ); // buffer protocol strides are in bytes, AMReX strides are elements auto const strides = py::make_tuple( // Strides (in bytes) for each index sizeof(T) * a4.jstride, sizeof(T) * a4.kstride, - sizeof(T)//, - //sizeof(T) * a4.nstride + sizeof(T), + sizeof(T) * a4.nstride ); d["data"] = py::make_tuple(long(a4.dataPtr()), false); d["typestr"] = py::format_descriptor::format(); @@ -117,7 +121,25 @@ void make_Array4(py::module &m, std::string typestr) ); }) */ + .def("contains", &Array4::contains) + //.def("__contains__", &Array4::contains) + + //.def("__setitem__", [](Array4 & a4, IntVect v, T const value){ a4(v[0], v[1], v[2]) = value; }) + .def("__setitem__", [](Array4 & a4, std::array key, T const value){ + a4(key[0], key[1], key[2], key[3]) = value; + }) + .def("__setitem__", [](Array4 & a4, std::array key, T const value){ + a4(key[0], key[1], key[2]) = value; + }) + + // __getitem__ ; + + // free standing C++ functions: + m.def("lbound", &lbound< Array4 >); + m.def("ubound", &ubound< Array4 >); + m.def("length", &length< Array4 >); + m.def("makePolymorphic", &makePolymorphic< Array4 >); } void init_Array4(py::module &m) { @@ -148,13 +170,4 @@ void init_Array4(py::module &m) { ) ; */ - - // free standing C++ functions: - /* - contains - lbound - ubound - length - makePolymorphic - */ } diff --git a/tests/test_array4.py b/tests/test_array4.py index 0c8f0a42..a83467d4 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -21,6 +21,7 @@ def test_array4(): # from numpy (also a non-owning view) x = np.ones((2, 3, 4,)) arr = amrex.Array4_double(x) + assert(arr.nComp == 1) x[1, 1, 1] = 42 # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable @@ -28,18 +29,16 @@ def test_array4(): # copy to numpy c_arr2np = np.array(arr, copy=True) - assert(c_arr2np.ndim == 3) - #assert(c_arr2np.ndim == 4) + assert(c_arr2np.ndim == 4) assert(c_arr2np.dtype == np.dtype("double")) - np.testing.assert_array_equal(x, c_arr2np) + np.testing.assert_array_equal(x, c_arr2np[:, :, :, 0]) assert(c_arr2np[1, 1, 1] == 42) # view to numpy v_arr2np = np.array(arr, copy=False) - assert(v_arr2np.ndim == 3) - #assert(c_arr2np.ndim == 4) + assert(c_arr2np.ndim == 4) assert(v_arr2np.dtype == np.dtype("double")) - np.testing.assert_array_equal(x, v_arr2np) + np.testing.assert_array_equal(x, v_arr2np[:, :, :, 0]) assert(v_arr2np[1, 1, 1] == 42) # change original buffer once more From f06fdcf477056aef220324d0a2804fdf62b48805 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 21:20:12 -0800 Subject: [PATCH 13/43] Iterate Box --- src/Base/Box.cpp | 73 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/src/Base/Box.cpp b/src/Base/Box.cpp index 022a3e62..37fd24a9 100644 --- a/src/Base/Box.cpp +++ b/src/Base/Box.cpp @@ -11,14 +11,64 @@ #include #include +#include namespace py = pybind11; using namespace amrex; +namespace +{ + /** A little Wrapper class to iterate an amrex::Box via + * amrex::Box::next(). + */ + struct Box3DConstIter { + Box m_box; + std::optional m_it; + + Box3DConstIter(Box const & bx) : m_box(bx) { + m_it = m_box.smallEnd(); + } + + Box3DConstIter& operator++() { + // from FABio_ascii::write + if (m_it <= m_box.bigEnd()) { + m_box.next(m_it.value()); + return *this; + } + else + { + m_it = std::nullopt; + return *this; + } + } + + bool operator==(Box3DConstIter const & other) const { + return other.m_it == m_it; + } + + Box3DConstIter begin() const + { + return Box3DConstIter(m_box); + } + Box3DConstIter end() const + { + auto it = Box3DConstIter(m_box); + it.m_it = std::nullopt; + return it; + } + + IntVect operator*() const + { + return m_it.value(); + } + }; +} void init_Box(py::module &m) { py::class_< Direction >(m, "Direction"); + + py::class_< Box >(m, "Box") .def("__repr__", [](Box const & b) { @@ -106,5 +156,28 @@ void init_Box(py::module &m) { // __getitem__ + /* iterate Box index space */ + .def("__iter__", + [](Box const & bx) { + auto box_iter = Box3DConstIter(bx); + return py::make_iterator(box_iter.begin(), box_iter.end()); + }, + // Essential: keep object alive while iterator exists + py::keep_alive<0, 1>() + ) + + .def("lbound", [](Box const &, Box const & other){ return lbound(other); }) + .def("ubound", [](Box const &, Box const & other){ return ubound(other); }) + .def("begin", [](Box const &, Box const & other){ return begin(other); }) + .def("end", [](Box const &, Box const & other){ return end(other); }) + // already an attribute + //.def("length", [](Box const &, Box const & other){ return length(other); }) ; + + // free standing C++ functions: + m.def("lbound", [](Box const & other){ return lbound(other); }); + m.def("ubound", [](Box const & other){ return ubound(other); }); + m.def("begin", [](Box const & other){ return begin(other); }); + m.def("end", [](Box const & other){ return end(other); }); + m.def("length", [](Box const & other){ return length(other); }); } From 285aaec5bc87cb7b8b740f1598f010d03a91430a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 21:20:23 -0800 Subject: [PATCH 14/43] Docs: Verbose Testing --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index 9bd94d02..90537b61 100644 --- a/README.md +++ b/README.md @@ -147,6 +147,9 @@ python3 -m pytest tests/test_intvect.py # Run a single test (useful during debugging) python3 -m pytest tests/test_intvect.py::test_iv_conversions + +# Check all print and be maximum verbose +python3 -m pytest -vvvv -s tests/test_multifab.py::test_mfab_simple ``` ### Build Options From 32c8b6bef7c44c504271d059e789b54a769cf307 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 21:20:41 -0800 Subject: [PATCH 15/43] FabArray: Iterate via Array4 --- src/Base/MultiFab.cpp | 11 +++++++++-- tests/test_multifab.py | 29 ++++++++++++++++++++++------- 2 files changed, 31 insertions(+), 9 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 3aefbd01..cdc289c9 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -28,13 +28,13 @@ void init_MultiFab(py::module &m) { .def("is_nodal", py::overload_cast< int >(&FabArrayBase::is_nodal, py::const_)) + .def_property_readonly("nComp", &FabArrayBase::nComp) .def_property_readonly("num_comp", &FabArrayBase::nComp) + .def_property_readonly("size", &FabArrayBase::size) ; py::class_< FArrayBox >(m, "FArrayBox"); - py::class_< FabArray, FabArrayBase >(m, "FabArray_FArrayBox"); - //py::class_< MFIterWrapper >(m, "MFIterWrapper"); py::class_< MFIter >(m, "MFIter") .def("__repr__", @@ -94,6 +94,13 @@ void init_MultiFab(py::module &m) { */ ; + py::class_< FabArray, FabArrayBase >(m, "FabArray_FArrayBox") + //.def("array", py::overload_cast< const MFIter& >(&FabArray::array)) + //.def("const_array", &FabArray::const_array) + .def("array", [](FabArray & fa, MFIter const & mfi) {return fa.array(mfi);}) + .def("const_array", [](FabArray & fa, MFIter const & mfi) {return fa.const_array(mfi);}) + ; + py::class_< MultiFab, FabArray >(m, "MultiFab") .def("__repr__", [](MultiFab const & mf) { diff --git a/tests/test_multifab.py b/tests/test_multifab.py index f0167788..2848ada9 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -10,13 +10,28 @@ def test_mfab_loop(mfab, nghost): for mfi in mfab: print(mfi) bx = mfi.tilebox() - #marr = mfab.array(mfi) # Array4: add __array_interface__ here - - #for i, j, k in bx: - # mar[i, j, k, 0] = 10.0 * i - # mar[i, j, k, 1] = 10.0 * j - # mar[i, j, k, 2] = 10.0 * k - # print(i,j,k) + marr = mfab.array(mfi) # Array4: add __array_interface__ here + + print(mfab) + print(mfab.num_comp) + print(mfab.size) + print(marr.size) + print(marr.nComp) + + # slow, index by index assignment + for i, j, k in bx: + #print(i,j,k) + marr[i, j, k] = 10.0 * i + #marr[i, j, k, 0] = 10.0 * i + #marr[i, j, k, 1] = 10.0 * j + #marr[i, j, k, 2] = 10.0 * k + + # fast, range based assignment + # TODO + + # challenge: offset from our index space for... + # extra test: numpy assignment + # extra test: cupy assignment def test_mfab_simple(mfab): From c58abd35e68ef5a937d28243db0af6d429e4dd92 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 22:44:16 -0800 Subject: [PATCH 16/43] Box: Fix Iteration Range, Add Shift/Grow --- src/Base/Box.cpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/Base/Box.cpp b/src/Base/Box.cpp index 37fd24a9..d53e6c87 100644 --- a/src/Base/Box.cpp +++ b/src/Base/Box.cpp @@ -4,6 +4,7 @@ * License: BSD-3-Clause-LBNL */ #include +#include #include #include @@ -31,7 +32,7 @@ namespace Box3DConstIter& operator++() { // from FABio_ascii::write - if (m_it <= m_box.bigEnd()) { + if (m_it < m_box.bigEnd()) { m_box.next(m_it.value()); return *this; } @@ -82,6 +83,8 @@ void init_Box(py::module &m) { .def(py::init< IntVect const &, IntVect const &, IntVect const & >()) //.def(py::init< IntVect const &, IntVect const &, IndexType >()) + .def_property_readonly("small_end", [](Box const & bx){ return bx.smallEnd(); }) + .def_property_readonly("big_end", [](Box const & bx){ return bx.bigEnd(); }) /* .def_property("small_end", &Box::smallEnd, @@ -89,7 +92,7 @@ void init_Box(py::module &m) { .def_property("big_end", &Box::bigEnd, &Box::setBig) - */ + .def_property("type", py::overload_cast<>(&Box::type, py::const_), &Box::setType) @@ -126,14 +129,22 @@ void init_Box(py::module &m) { // atOffset // atOffset3d // setRange - // shift // shiftHalf + */ + .def("shift", [](Box & bx, IntVect const& iv) { return bx.shift(iv); }) + + .def(py::self + IntVect()) + .def(py::self - IntVect()) + .def(py::self += IntVect()) + .def(py::self -= IntVect()) .def("convert", py::overload_cast< IndexType >(&Box::convert)) .def("convert", py::overload_cast< IntVect const & >(&Box::convert)) + .def("grow", [](Box & bx, IntVect const& iv) { return bx.grow(iv); }) + //.def("surrounding_nodes", // py::overload_cast< >(&Box::surroundingNodes)) //.def("surrounding_nodes", From 3a3590a9d3559961768d43fd780eea661abeb25d Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 22:44:50 -0800 Subject: [PATCH 17/43] Array4: Setter/Getter --- src/Base/Array4.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index da8e8f72..d5728143 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -124,15 +124,22 @@ void make_Array4(py::module &m, std::string typestr) .def("contains", &Array4::contains) //.def("__contains__", &Array4::contains) - //.def("__setitem__", [](Array4 & a4, IntVect v, T const value){ a4(v[0], v[1], v[2]) = value; }) - .def("__setitem__", [](Array4 & a4, std::array key, T const value){ + // setter & getter + .def("__setitem__", [](Array4 & a4, IntVect const & v, T const value){ a4(v) = value; }) + .def("__setitem__", [](Array4 & a4, std::array const key, T const value){ a4(key[0], key[1], key[2], key[3]) = value; }) - .def("__setitem__", [](Array4 & a4, std::array key, T const value){ + .def("__setitem__", [](Array4 & a4, std::array const key, T const value){ a4(key[0], key[1], key[2]) = value; }) - // __getitem__ + .def("__getitem__", [](Array4 & a4, IntVect const & v){ return a4(v); }) + .def("__getitem__", [](Array4 & a4, std::array const key){ + return a4(key[0], key[1], key[2], key[3]); + }) + .def("__getitem__", [](Array4 & a4, std::array const key){ + return a4(key[0], key[1], key[2]); + }) ; // free standing C++ functions: From 48612b7b379c0fd19b88c9686154da94b2fe6b48 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 8 Mar 2022 22:45:12 -0800 Subject: [PATCH 18/43] MultiFab: Numpy Roundtrip --- src/Base/MultiFab.cpp | 2 ++ tests/test_multifab.py | 55 +++++++++++++++++++++++++++++------------- 2 files changed, 40 insertions(+), 17 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index cdc289c9..35db1e4c 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -31,6 +31,8 @@ void init_MultiFab(py::module &m) { .def_property_readonly("nComp", &FabArrayBase::nComp) .def_property_readonly("num_comp", &FabArrayBase::nComp) .def_property_readonly("size", &FabArrayBase::size) + + .def_property_readonly("nGrowVect", &FabArrayBase::nGrowVect) ; py::class_< FArrayBox >(m, "FArrayBox"); diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 2848ada9..48b0f7df 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -7,31 +7,52 @@ @pytest.mark.parametrize("nghost", [0, 1]) def test_mfab_loop(mfab, nghost): + ngv = mfab.nGrowVect + print(ngv) + for mfi in mfab: print(mfi) - bx = mfi.tilebox() - marr = mfab.array(mfi) # Array4: add __array_interface__ here + bx = mfi.tilebox().grow(ngv) + marr = mfab.array(mfi) - print(mfab) - print(mfab.num_comp) - print(mfab.size) - print(marr.size) - print(marr.nComp) + #print(mfab) + #print(mfab.num_comp) + #print(mfab.size) + #print(marr.size) + #print(marr.nComp) # slow, index by index assignment - for i, j, k in bx: - #print(i,j,k) - marr[i, j, k] = 10.0 * i - #marr[i, j, k, 0] = 10.0 * i - #marr[i, j, k, 1] = 10.0 * j - #marr[i, j, k, 2] = 10.0 * k + three_comps = mfab.num_comp == 3 + if three_comps: + for i, j, k in bx: + #print(i,j,k) + marr[i, j, k, 0] = 10.0 * i + marr[i, j, k, 1] = 10.0 * j + marr[i, j, k, 2] = 10.0 * k + else: + for i, j, k in bx: + #print(i,j,k) + marr[i, j, k] = 10.0 * i # fast, range based assignment - # TODO + # challenge: offset from index space + #bx_zeroshift = bx - bx.small_end - mfab.nGrowVect + + # numpy assignment: including guard/ghost region + marr_np = np.array(marr, copy=False) + print(marr_np.shape) + + assert(marr_np[0, 0, 0, 0] == marr[bx.small_end]) + # assert(marr_np[-1, -1, -1, -1] == marr[bx.big_end]) # FIXME + + #marr_np[24:200, :, :, :] = 42. # this should fail + #marr_np[:, :, :] = 42. + marr_np[:, :, :, :] = 42. + assert(marr_np[0, 0, 0, 0] == marr[bx.small_end]) + assert(marr_np[-1, -1, -1, -1] == marr[bx.big_end]) - # challenge: offset from our index space for... - # extra test: numpy assignment - # extra test: cupy assignment + # separate test: cupy assignment & reading + # TODO def test_mfab_simple(mfab): From 47cdbb11f2e70c97e5a99d5eddedf0f90d51e170 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 10 Mar 2022 00:09:30 -0800 Subject: [PATCH 19/43] fix: test_multifab is_nodal only positive UBSAN: ``` tests/test_multifab.py::test_mfab_simple[mfab0] /home/axel/src/pyamrex/build/_deps/fetchedamrex-src/Src/Base/AMReX_IndexType.H:136:48: runtime error: shift exponent -1 is negative ``` --- tests/test_multifab.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 48b0f7df..4df12a2d 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -57,7 +57,8 @@ def test_mfab_loop(mfab, nghost): def test_mfab_simple(mfab): assert(mfab.is_all_cell_centered) - assert(all(not mfab.is_nodal(i) for i in [-1, 0, 1, 2])) + #assert(all(not mfab.is_nodal(i) for i in [-1, 0, 1, 2])) # -1?? + assert(all(not mfab.is_nodal(i) for i in [0, 1, 2])) for i in range(mfab.num_comp): mfab.set_val(-10 * (i + 1), i, 1) From c7e92e422b6592f5d2d064b0f80b780a0bd18f43 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 10 Mar 2022 00:11:36 -0800 Subject: [PATCH 20/43] Win CI: More Verbose --- .github/workflows/windows.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 5bdce2aa..c2b27190 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -19,7 +19,7 @@ jobs: - name: Unit tests shell: cmd run: | - python3 -m pytest tests + python3 -m pytest -s -vvvv tests clang: name: Clang w/o MPI From 8e228420cb3f96539d9e1418dbe421a53b2435f6 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 24 Mar 2022 12:24:28 -0700 Subject: [PATCH 21/43] Comment: Segfault on Windows: Numpy to Array4 --- tests/test_array4.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_array4.py b/tests/test_array4.py index a83467d4..08366ab0 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -28,7 +28,7 @@ def test_array4(): # assert(arr[1, 1, 1] == 42) # copy to numpy - c_arr2np = np.array(arr, copy=True) + c_arr2np = np.array(arr, copy=True) # segfaults on Windows assert(c_arr2np.ndim == 4) assert(c_arr2np.dtype == np.dtype("double")) np.testing.assert_array_equal(x, c_arr2np[:, :, :, 0]) From 3f8ec99b829bc0f3005395b9bf9b5426baf35d4e Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 24 Mar 2022 14:56:16 -0700 Subject: [PATCH 22/43] Print: __array_interface__ in test --- tests/test_array4.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_array4.py b/tests/test_array4.py index 08366ab0..768c5b44 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -20,7 +20,9 @@ def test_array4_empty(): def test_array4(): # from numpy (also a non-owning view) x = np.ones((2, 3, 4,)) + print(x.__array_interface__) arr = amrex.Array4_double(x) + print(arr.__array_interface__) assert(arr.nComp == 1) x[1, 1, 1] = 42 From 8d17e992a5895bf7455c5a0da938700087f0a673 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 09:06:28 -0700 Subject: [PATCH 23/43] Windows: Verbose Array4 Tests --- .github/workflows/windows.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index c2b27190..08ac5f32 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -19,7 +19,7 @@ jobs: - name: Unit tests shell: cmd run: | - python3 -m pytest -s -vvvv tests + python3 -m pytest -s -vvvv tests/test_array4.py clang: name: Clang w/o MPI @@ -42,7 +42,7 @@ jobs: python3 -m pip install -v . - name: Unit tests run: | - ctest --test-dir build -C Release --output-on-failure + python3 -m pytest -s -vvvv tests/test_array4.py # C:\\hostedtoolcache\\windows\\python\\3.9.10\\x64\\python3.exe # C:\\hostedtoolcache\\windows\\Python\\3.10.2\\x64\\python3.exe From d9ee624fa21299ffc619e70eb11187203b71ca45 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 09:30:08 -0700 Subject: [PATCH 24/43] `pytest -s` not working on Win? --- tests/test_array4.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/test_array4.py b/tests/test_array4.py index 768c5b44..993ddcea 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -17,12 +17,14 @@ def test_array4_empty(): assert(emptyc.size == 0) assert(emptyc.nComp == 0) -def test_array4(): +def test_array4(capsys): # from numpy (also a non-owning view) x = np.ones((2, 3, 4,)) - print(x.__array_interface__) + with capsys.disabled(): + print(x.__array_interface__) arr = amrex.Array4_double(x) - print(arr.__array_interface__) + with capsys.disabled(): + print(arr.__array_interface__) assert(arr.nComp == 1) x[1, 1, 1] = 42 From e83e8eb53d2a3a999a9ec81fd6ec22a146cf83ac Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 09:43:17 -0700 Subject: [PATCH 25/43] Return before access violation --- tests/test_array4.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_array4.py b/tests/test_array4.py index 993ddcea..f65dd140 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -26,6 +26,7 @@ def test_array4(capsys): with capsys.disabled(): print(arr.__array_interface__) assert(arr.nComp == 1) + return x[1, 1, 1] = 42 # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable From d7290876f2d7051584266af6dc6fcae17f34c518 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 15:03:05 -0700 Subject: [PATCH 26/43] Update MultiFab Test: Component Access --- tests/test_multifab.py | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index 4df12a2d..d702f36e 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -8,10 +8,9 @@ @pytest.mark.parametrize("nghost", [0, 1]) def test_mfab_loop(mfab, nghost): ngv = mfab.nGrowVect - print(ngv) + print(f"\n mfab={mfab}, mfab.nGrowVect={ngv}") for mfi in mfab: - print(mfi) bx = mfi.tilebox().grow(ngv) marr = mfab.array(mfi) @@ -21,7 +20,11 @@ def test_mfab_loop(mfab, nghost): #print(marr.size) #print(marr.nComp) - # slow, index by index assignment + # index by index assignment + # notes: + # - this is AMReX Array4, F-order indices + # - even though we iterate by fastest varying index, + # such loops are naturally very slow in Python three_comps = mfab.num_comp == 3 if three_comps: for i, j, k in bx: @@ -34,19 +37,30 @@ def test_mfab_loop(mfab, nghost): #print(i,j,k) marr[i, j, k] = 10.0 * i - # fast, range based assignment - # challenge: offset from index space - #bx_zeroshift = bx - bx.small_end - mfab.nGrowVect + # note: offset from index space in numpy + # in numpy, we start indices from zero, not small_end - # numpy assignment: including guard/ghost region + # numpy representation: non-copying view, including the + # guard/ghost region + # note: in numpy, indices are in C-order! marr_np = np.array(marr, copy=False) - print(marr_np.shape) + # check the values at start/end are the same: first component assert(marr_np[0, 0, 0, 0] == marr[bx.small_end]) - # assert(marr_np[-1, -1, -1, -1] == marr[bx.big_end]) # FIXME - - #marr_np[24:200, :, :, :] = 42. # this should fail - #marr_np[:, :, :] = 42. + assert(marr_np[0, -1, -1, -1] == marr[bx.big_end]) + # same check, but for all components + for n in range(mfab.num_comp): + small_end_comp = list(bx.small_end) + [n] + big_end_comp = list(bx.big_end) + [n] + assert(marr_np[n, 0, 0, 0] == marr[small_end_comp]) + assert(marr_np[n, -1, -1, -1] == marr[big_end_comp]) + + # now we do some faster assignments, using range based access + # this should fail as out-of-bounds, but does not + # does Numpy not check array access for non-owned views? + #marr_np[24:200, :, :, :] = 42. + + # all components and all indices set at once to 42 marr_np[:, :, :, :] = 42. assert(marr_np[0, 0, 0, 0] == marr[bx.small_end]) assert(marr_np[-1, -1, -1, -1] == marr[bx.big_end]) From 56dfeb36fafec280b7ec0a56dfcc1b3b94fcc5e1 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 11:35:13 -0700 Subject: [PATCH 27/43] Fix: Array4 <-> numpy/buffer protocol --- src/Base/Array4.cpp | 37 +++++++++++++++++++++++-------------- tests/test_array4.py | 15 ++++++++------- 2 files changed, 31 insertions(+), 21 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index d5728143..a13661b1 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -7,8 +7,8 @@ #include #include -#include #include +#include #include #include @@ -50,41 +50,50 @@ void make_Array4(py::module &m, std::string typestr) .def(py::init([](py::array_t & arr) { py::buffer_info buf = arr.request(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(buf.ndim == 3, + "We can only create amrex::Array4 views into 3D Python arrays at the moment."); + auto a4 = std::make_unique< Array4 >(); a4.get()->p = (T*)buf.ptr; a4.get()->begin = Dim3{0, 0, 0}; - // TODO: likely C->F index conversion here + // C->F index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - a4.get()->end.x = (int)buf.shape.at(0); + a4.get()->end.x = (int)buf.shape.at(2); a4.get()->end.y = (int)buf.shape.at(1); - a4.get()->end.z = (int)buf.shape.at(2); + a4.get()->end.z = (int)buf.shape.at(0); a4.get()->ncomp = 1; // buffer protocol strides are in bytes, AMReX strides are elements - a4.get()->jstride = (int)buf.strides.at(0) / sizeof(T); - a4.get()->kstride = (int)buf.strides.at(1) / sizeof(T); - a4.get()->nstride = (int)buf.strides.at(2) * (int)buf.shape.at(2) / sizeof(T); + a4.get()->jstride = (int)buf.strides.at(1) / sizeof(T); + a4.get()->kstride = (int)buf.strides.at(0) / sizeof(T); + a4.get()->nstride = (int)buf.strides.at(0) * (int)buf.shape.at(0) / sizeof(T); + + + std::cout << "(int)buf.strides.at(0)=" << (int)buf.strides.at(0) << std::endl; + std::cout << "(int)buf.strides.at(1)=" << (int)buf.strides.at(1) << std::endl; + std::cout << "(int)buf.strides.at(2)=" << (int)buf.strides.at(2) << std::endl; + return a4; })) .def_property_readonly("__array_interface__", [](Array4 const & a4) { auto d = py::dict(); auto const len = length(a4); - // TODO: likely F->C index conversion here + // F->C index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; py::print("ncomp"); py::print(a4.ncomp); auto shape = py::make_tuple( // Buffer dimensions - len.x < 0 ? 0 : len.x, + a4.ncomp, + len.z < 0 ? 0 : len.z, len.y < 0 ? 0 : len.y, - len.z < 0 ? 0 : len.z, // zero-size shall not have negative dimension - a4.ncomp + len.x < 0 ? 0 : len.x // zero-size shall not have negative dimension ); // buffer protocol strides are in bytes, AMReX strides are elements auto const strides = py::make_tuple( // Strides (in bytes) for each index - sizeof(T) * a4.jstride, + sizeof(T) * a4.nstride, sizeof(T) * a4.kstride, - sizeof(T), - sizeof(T) * a4.nstride + sizeof(T) * a4.jstride, + sizeof(T) ); d["data"] = py::make_tuple(long(a4.dataPtr()), false); d["typestr"] = py::format_descriptor::format(); diff --git a/tests/test_array4.py b/tests/test_array4.py index f65dd140..90e72602 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -26,7 +26,6 @@ def test_array4(capsys): with capsys.disabled(): print(arr.__array_interface__) assert(arr.nComp == 1) - return x[1, 1, 1] = 42 # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable @@ -36,25 +35,27 @@ def test_array4(capsys): c_arr2np = np.array(arr, copy=True) # segfaults on Windows assert(c_arr2np.ndim == 4) assert(c_arr2np.dtype == np.dtype("double")) - np.testing.assert_array_equal(x, c_arr2np[:, :, :, 0]) - assert(c_arr2np[1, 1, 1] == 42) + with capsys.disabled(): + print(c_arr2np.__array_interface__) + np.testing.assert_array_equal(x, c_arr2np[0, :, :, :]) + assert(c_arr2np[0, 1, 1, 1] == 42) # view to numpy v_arr2np = np.array(arr, copy=False) assert(c_arr2np.ndim == 4) assert(v_arr2np.dtype == np.dtype("double")) - np.testing.assert_array_equal(x, v_arr2np[:, :, :, 0]) - assert(v_arr2np[1, 1, 1] == 42) + np.testing.assert_array_equal(x, v_arr2np[0, :, :, :]) + assert(v_arr2np[0, 1, 1, 1] == 42) # change original buffer once more x[1, 1, 1] = 43 - assert(v_arr2np[1, 1, 1] == 43) + assert(v_arr2np[0, 1, 1, 1] == 43) # copy array4 (view) c_arr = amrex.Array4_double(arr) v_carr2np = np.array(c_arr, copy=False) x[1, 1, 1] = 44 - assert(v_carr2np[1, 1, 1] == 44) + assert(v_carr2np[0, 1, 1, 1] == 44) # from cupy From de77affdffe7152167976f50cd3d6261b8c0fbcc Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 14:09:40 -0700 Subject: [PATCH 28/43] Array4: 2D/1D note --- src/Base/Array4.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index a13661b1..c054a65d 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -52,6 +52,9 @@ void make_Array4(py::module &m, std::string typestr) AMREX_ALWAYS_ASSERT_WITH_MESSAGE(buf.ndim == 3, "We can only create amrex::Array4 views into 3D Python arrays at the moment."); + // TODO: + // In 2D, Array4 still needs to be accessed with (i,j,k) or (i,j,k,n), with k = 0. + // Likewise in 1D. auto a4 = std::make_unique< Array4 >(); a4.get()->p = (T*)buf.ptr; From 8ba5a0aa9a523c80231d719f2e4a5891e75ef243 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 15:01:49 -0700 Subject: [PATCH 29/43] Array4: 3D set nstride 3D == no component: stride should not matter --- src/Base/Array4.cpp | 45 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index c054a65d..38f7af35 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -55,26 +55,31 @@ void make_Array4(py::module &m, std::string typestr) // TODO: // In 2D, Array4 still needs to be accessed with (i,j,k) or (i,j,k,n), with k = 0. // Likewise in 1D. + // We could also add support for 4D numpy arrays, treating the slowest + // varying index as component "n". auto a4 = std::make_unique< Array4 >(); a4.get()->p = (T*)buf.ptr; a4.get()->begin = Dim3{0, 0, 0}; // C->F index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - a4.get()->end.x = (int)buf.shape.at(2); + a4.get()->end.x = (int)buf.shape.at(2); // fastest varying index a4.get()->end.y = (int)buf.shape.at(1); a4.get()->end.z = (int)buf.shape.at(0); a4.get()->ncomp = 1; // buffer protocol strides are in bytes, AMReX strides are elements - a4.get()->jstride = (int)buf.strides.at(1) / sizeof(T); + a4.get()->jstride = (int)buf.strides.at(1) / sizeof(T); // fastest varying index a4.get()->kstride = (int)buf.strides.at(0) / sizeof(T); - a4.get()->nstride = (int)buf.strides.at(0) * (int)buf.shape.at(0) / sizeof(T); + // 3D == no component: stride here should not matter + a4.get()->nstride = a4.get()->kstride * (int)buf.shape.at(0); std::cout << "(int)buf.strides.at(0)=" << (int)buf.strides.at(0) << std::endl; std::cout << "(int)buf.strides.at(1)=" << (int)buf.strides.at(1) << std::endl; std::cout << "(int)buf.strides.at(2)=" << (int)buf.strides.at(2) << std::endl; + // todo: we could check and store here if the array buffer we got is read-only + return a4; })) @@ -83,30 +88,50 @@ void make_Array4(py::module &m, std::string typestr) auto const len = length(a4); // F->C index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - py::print("ncomp"); - py::print(a4.ncomp); - auto shape = py::make_tuple( // Buffer dimensions + //py::print("ncomp"); + //py::print(a4.ncomp); + // Buffer dimensions: zero-size shall not have negative dimension + auto shape = py::make_tuple( a4.ncomp, len.z < 0 ? 0 : len.z, len.y < 0 ? 0 : len.y, - len.x < 0 ? 0 : len.x // zero-size shall not have negative dimension + len.x < 0 ? 0 : len.x // fastest varying index ); // buffer protocol strides are in bytes, AMReX strides are elements auto const strides = py::make_tuple( // Strides (in bytes) for each index sizeof(T) * a4.nstride, sizeof(T) * a4.kstride, sizeof(T) * a4.jstride, - sizeof(T) + sizeof(T) // fastest varying index ); - d["data"] = py::make_tuple(long(a4.dataPtr()), false); - d["typestr"] = py::format_descriptor::format(); + bool const read_only = false; + d["data"] = py::make_tuple(long(a4.dataPtr()), read_only); + //d["offset"] = 0; + //d["mask"] = py::none(); + d["shape"] = shape; d["strides"] = strides; + // we could set this after checking the strides are C-style contiguous: // d["strides"] = py::none(); + + d["typestr"] = py::format_descriptor::format(); d["version"] = 3; return d; }) + + // TODO: __cuda_array_interface__ + // https://numba.readthedocs.io/en/latest/cuda/cuda_array_interface.html + + + // TODO: __dlpack__ + // DLPack protocol (CPU, NVIDIA GPU, AMD GPU, Intel GPU, etc.) + // https://dmlc.github.io/dlpack/latest/ + // https://data-apis.org/array-api/latest/design_topics/data_interchange.html + // https://github.com/data-apis/consortium-feedback/issues/1 + // https://github.com/dmlc/dlpack/blob/master/include/dlpack/dlpack.h + + // not sure if useful to have this implemented on top /* .def_buffer([](Array4 & a4) -> py::buffer_info { From 9db3214aa4c3c33047752ef5bdb3d88ebcc3e747 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 15:25:33 -0700 Subject: [PATCH 30/43] hack: early return --- tests/test_array4.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/test_array4.py b/tests/test_array4.py index 90e72602..737db7d3 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -22,6 +22,7 @@ def test_array4(capsys): x = np.ones((2, 3, 4,)) with capsys.disabled(): print(x.__array_interface__) + print(x.dtype) arr = amrex.Array4_double(x) with capsys.disabled(): print(arr.__array_interface__) @@ -30,6 +31,9 @@ def test_array4(capsys): x[1, 1, 1] = 42 # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable # assert(arr[1, 1, 1] == 42) + + # hack: + return # copy to numpy c_arr2np = np.array(arr, copy=True) # segfaults on Windows From ae07d7c9164f882ed6335d07fa05e37f77a73bc4 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 16:19:57 -0700 Subject: [PATCH 31/43] do not return --- tests/test_array4.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/tests/test_array4.py b/tests/test_array4.py index 737db7d3..76941b51 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -31,9 +31,6 @@ def test_array4(capsys): x[1, 1, 1] = 42 # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable # assert(arr[1, 1, 1] == 42) - - # hack: - return # copy to numpy c_arr2np = np.array(arr, copy=True) # segfaults on Windows From 9262e7dac93e9987599123656000f9059d4a6611 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 22:31:16 -0700 Subject: [PATCH 32/43] Array4: Check Buffer Type --- src/Base/Array4.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index 38f7af35..c19ca87e 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -58,8 +58,13 @@ void make_Array4(py::module &m, std::string typestr) // We could also add support for 4D numpy arrays, treating the slowest // varying index as component "n". + if (buf.format != py::format_descriptor::format()) + throw std::runtime_error("Incompatible format: expected '" + + py::format_descriptor::format() + + "' and received '" + buf.format + "'!"); + auto a4 = std::make_unique< Array4 >(); - a4.get()->p = (T*)buf.ptr; + a4.get()->p = static_cast(buf.ptr); a4.get()->begin = Dim3{0, 0, 0}; // C->F index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; From 699345655bd51aa73675bbff7e78b362eddf680c Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 22:36:50 -0700 Subject: [PATCH 33/43] Array4: remove py::buffer_protocol --- src/Base/Array4.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index c19ca87e..948d6fdd 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -23,7 +23,7 @@ void make_Array4(py::module &m, std::string typestr) { // dispatch simpler via: py::format_descriptor::format() naming auto const array_name = std::string("Array4_").append(typestr); - py::class_< Array4 >(m, array_name.c_str(), py::buffer_protocol()) + py::class_< Array4 >(m, array_name.c_str() /*, py::buffer_protocol() */) .def("__repr__", [](Array4 const & a4) { std::stringstream s; From 33847e7da7780061cf18696b37e56d3445ed423f Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 23:20:01 -0700 Subject: [PATCH 34/43] __array_interface__ data as std::intptr_t --- src/Base/Array4.cpp | 3 ++- tests/test_array4.py | 9 +++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index 948d6fdd..d75e74b7 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include @@ -110,7 +111,7 @@ void make_Array4(py::module &m, std::string typestr) sizeof(T) // fastest varying index ); bool const read_only = false; - d["data"] = py::make_tuple(long(a4.dataPtr()), read_only); + d["data"] = py::make_tuple(std::intptr_t(a4.dataPtr()), read_only); //d["offset"] = 0; //d["mask"] = py::none(); diff --git a/tests/test_array4.py b/tests/test_array4.py index 76941b51..0ee45a2a 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -28,9 +28,14 @@ def test_array4(capsys): print(arr.__array_interface__) assert(arr.nComp == 1) + # change original array x[1, 1, 1] = 42 - # TypeError: 'amrex.amrex_pybind.Array4_double' object is not subscriptable - # assert(arr[1, 1, 1] == 42) + # check values in Array4 view changed + assert(arr[1, 1, 1] == 42) + assert(arr[1, 1, 1, 0] == 42) # with component + # check existing values stayed + assert(arr[0, 0, 0] == 1) + assert(arr[3, 2, 1] == 1) # copy to numpy c_arr2np = np.array(arr, copy=True) # segfaults on Windows From b803492c90e99278a5b2cb2aa5e600080c1c9f25 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 23:50:36 -0700 Subject: [PATCH 35/43] Array4: Cleanup --- src/Base/Array4.cpp | 19 +++++++------------ tests/test_array4.py | 12 ++++-------- 2 files changed, 11 insertions(+), 20 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index d75e74b7..7e181528 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -79,11 +79,6 @@ void make_Array4(py::module &m, std::string typestr) // 3D == no component: stride here should not matter a4.get()->nstride = a4.get()->kstride * (int)buf.shape.at(0); - - std::cout << "(int)buf.strides.at(0)=" << (int)buf.strides.at(0) << std::endl; - std::cout << "(int)buf.strides.at(1)=" << (int)buf.strides.at(1) << std::endl; - std::cout << "(int)buf.strides.at(2)=" << (int)buf.strides.at(2) << std::endl; - // todo: we could check and store here if the array buffer we got is read-only return a4; @@ -94,8 +89,6 @@ void make_Array4(py::module &m, std::string typestr) auto const len = length(a4); // F->C index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - //py::print("ncomp"); - //py::print(a4.ncomp); // Buffer dimensions: zero-size shall not have negative dimension auto shape = py::make_tuple( a4.ncomp, @@ -104,7 +97,7 @@ void make_Array4(py::module &m, std::string typestr) len.x < 0 ? 0 : len.x // fastest varying index ); // buffer protocol strides are in bytes, AMReX strides are elements - auto const strides = py::make_tuple( // Strides (in bytes) for each index + auto const strides = py::make_tuple( sizeof(T) * a4.nstride, sizeof(T) * a4.kstride, sizeof(T) * a4.jstride, @@ -112,13 +105,15 @@ void make_Array4(py::module &m, std::string typestr) ); bool const read_only = false; d["data"] = py::make_tuple(std::intptr_t(a4.dataPtr()), read_only); - //d["offset"] = 0; - //d["mask"] = py::none(); + //d["offset"] = 0; // default + //d["mask"] = py::none(); // default d["shape"] = shape; + // we could also set this after checking the strides are C-style contiguous: + //if (is_contiguous(shape, strides)) + // d["strides"] = py::none(); // C-style contiguous + //else d["strides"] = strides; - // we could set this after checking the strides are C-style contiguous: - // d["strides"] = py::none(); d["typestr"] = py::format_descriptor::format(); d["version"] = 3; diff --git a/tests/test_array4.py b/tests/test_array4.py index 0ee45a2a..7a94b770 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -17,15 +17,12 @@ def test_array4_empty(): assert(emptyc.size == 0) assert(emptyc.nComp == 0) -def test_array4(capsys): +def test_array4(): # from numpy (also a non-owning view) x = np.ones((2, 3, 4,)) - with capsys.disabled(): - print(x.__array_interface__) - print(x.dtype) + print(f"\nx: {x.__array_interface__} {x.dtype}") arr = amrex.Array4_double(x) - with capsys.disabled(): - print(arr.__array_interface__) + print(f"arr: {arr.__array_interface__}") assert(arr.nComp == 1) # change original array @@ -41,8 +38,7 @@ def test_array4(capsys): c_arr2np = np.array(arr, copy=True) # segfaults on Windows assert(c_arr2np.ndim == 4) assert(c_arr2np.dtype == np.dtype("double")) - with capsys.disabled(): - print(c_arr2np.__array_interface__) + print(f"c_arr2np: {c_arr2np.__array_interface__}") np.testing.assert_array_equal(x, c_arr2np[0, :, :, :]) assert(c_arr2np[0, 1, 1, 1] == 42) From 9615236bac8454eb430a066291c1c39d9405aaa2 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Fri, 25 Mar 2022 23:50:55 -0700 Subject: [PATCH 36/43] Windows CI: Cleanup --- .github/workflows/windows.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 08ac5f32..5bdce2aa 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -19,7 +19,7 @@ jobs: - name: Unit tests shell: cmd run: | - python3 -m pytest -s -vvvv tests/test_array4.py + python3 -m pytest tests clang: name: Clang w/o MPI @@ -42,7 +42,7 @@ jobs: python3 -m pip install -v . - name: Unit tests run: | - python3 -m pytest -s -vvvv tests/test_array4.py + ctest --test-dir build -C Release --output-on-failure # C:\\hostedtoolcache\\windows\\python\\3.9.10\\x64\\python3.exe # C:\\hostedtoolcache\\windows\\Python\\3.10.2\\x64\\python3.exe From 342b1d2171affb9c6b088338445f472f1a55c530 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 00:16:41 -0700 Subject: [PATCH 37/43] Array4: Always at least size 1 per dim --- src/Base/Array4.cpp | 39 ++++++++------------------------------- 1 file changed, 8 insertions(+), 31 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index 7e181528..4e3f1f14 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -24,7 +24,7 @@ void make_Array4(py::module &m, std::string typestr) { // dispatch simpler via: py::format_descriptor::format() naming auto const array_name = std::string("Array4_").append(typestr); - py::class_< Array4 >(m, array_name.c_str() /*, py::buffer_protocol() */) + py::class_< Array4 >(m, array_name.c_str()) .def("__repr__", [](Array4 const & a4) { std::stringstream s; @@ -89,12 +89,12 @@ void make_Array4(py::module &m, std::string typestr) auto const len = length(a4); // F->C index conversion here // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - // Buffer dimensions: zero-size shall not have negative dimension + // Buffer dimensions: zero-size shall not skip dimension auto shape = py::make_tuple( a4.ncomp, - len.z < 0 ? 0 : len.z, - len.y < 0 ? 0 : len.y, - len.x < 0 ? 0 : len.x // fastest varying index + len.z <= 0 ? 1 : len.z, + len.y <= 0 ? 1 : len.y, + len.x <= 0 ? 1 : len.x // fastest varying index ); // buffer protocol strides are in bytes, AMReX strides are elements auto const strides = py::make_tuple( @@ -105,6 +105,9 @@ void make_Array4(py::module &m, std::string typestr) ); bool const read_only = false; d["data"] = py::make_tuple(std::intptr_t(a4.dataPtr()), read_only); + // note: if we want to keep the same global indexing with non-zero + // box small_end as in AMReX, then we can explore playing with + // this offset as well //d["offset"] = 0; // default //d["mask"] = py::none(); // default @@ -133,32 +136,6 @@ void make_Array4(py::module &m, std::string typestr) // https://github.com/dmlc/dlpack/blob/master/include/dlpack/dlpack.h - // not sure if useful to have this implemented on top -/* - .def_buffer([](Array4 & a4) -> py::buffer_info { - auto const len = length(a4); - // TODO: likely F->C index conversion here - // p[(i-begin.x)+(j-begin.y)*jstride+(k-begin.z)*kstride+n*nstride]; - auto shape = { // Buffer dimensions - len.x < 0 ? 0 : len.x, - len.y < 0 ? 0 : len.y, - len.z < 0 ? 0 : len.z//, // zero-size shall not have negative dimension - //a4.ncomp - }; - // buffer protocol strides are in bytes, AMReX strides are elements - auto const strides = { // Strides (in bytes) for each index - sizeof(T) * a4.jstride, - sizeof(T) * a4.kstride, - sizeof(T)//, - //sizeof(T) * a4.nstride - }; - return py::buffer_info( - a4.dataPtr(), - shape, - strides - ); - }) -*/ .def("contains", &Array4::contains) //.def("__contains__", &Array4::contains) From 19b9e786fabb5cbb95e643425e9674f4e583a4c9 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 00:35:04 -0700 Subject: [PATCH 38/43] README: No-Capture & Verbose Reword --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 90537b61..a13bbb4a 100644 --- a/README.md +++ b/README.md @@ -148,8 +148,8 @@ python3 -m pytest tests/test_intvect.py # Run a single test (useful during debugging) python3 -m pytest tests/test_intvect.py::test_iv_conversions -# Check all print and be maximum verbose -python3 -m pytest -vvvv -s tests/test_multifab.py::test_mfab_simple +# Run all tests, do not capture "print" output and be verbose +python3 -m pytest -s -vvvv tests/ ``` ### Build Options From a62e8f01d2b3cf4d48f10ae91af669c701773059 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 00:38:32 -0700 Subject: [PATCH 39/43] Update: Year 2021-2022 --- src/Base/Array4.cpp | 2 +- src/Base/Box.cpp | 2 +- src/Base/BoxArray.cpp | 2 +- src/Base/DistributionMapping.cpp | 2 +- src/Base/MultiFab.cpp | 4 ++-- src/pyAMReX.cpp | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index 4e3f1f14..d460eb5e 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * * Authors: Axel Huebl * License: BSD-3-Clause-LBNL diff --git a/src/Base/Box.cpp b/src/Base/Box.cpp index d53e6c87..6dbbcab2 100644 --- a/src/Base/Box.cpp +++ b/src/Base/Box.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * * Authors: Axel Huebl * License: BSD-3-Clause-LBNL diff --git a/src/Base/BoxArray.cpp b/src/Base/BoxArray.cpp index e63ab6f6..802b6da5 100644 --- a/src/Base/BoxArray.cpp +++ b/src/Base/BoxArray.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * * Authors: Axel Huebl * License: BSD-3-Clause-LBNL diff --git a/src/Base/DistributionMapping.cpp b/src/Base/DistributionMapping.cpp index 5f57d0cc..2a32a55d 100644 --- a/src/Base/DistributionMapping.cpp +++ b/src/Base/DistributionMapping.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * * Authors: Axel Huebl * License: BSD-3-Clause-LBNL diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index 35db1e4c..aa311bf6 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -1,6 +1,6 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * - * Authors: Axel Huebl, ... + * Authors: Axel Huebl * License: BSD-3-Clause-LBNL */ #include diff --git a/src/pyAMReX.cpp b/src/pyAMReX.cpp index 9a7987c5..22911b59 100644 --- a/src/pyAMReX.cpp +++ b/src/pyAMReX.cpp @@ -1,4 +1,4 @@ -/* Copyright 2021 The AMReX Community +/* Copyright 2021-2022 The AMReX Community * * Authors: Axel Huebl * License: BSD-3-Clause-LBNL From 825bc2deeaa90b99232b72d1d461bcaf25087eff Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 00:40:38 -0700 Subject: [PATCH 40/43] MultiFab: Remove Old Iterator Draft Resolved more gentle already, just dead code left. --- src/Base/MultiFab.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/Base/MultiFab.cpp b/src/Base/MultiFab.cpp index aa311bf6..be4554f5 100644 --- a/src/Base/MultiFab.cpp +++ b/src/Base/MultiFab.cpp @@ -37,7 +37,6 @@ void init_MultiFab(py::module &m) { py::class_< FArrayBox >(m, "FArrayBox"); - //py::class_< MFIterWrapper >(m, "MFIterWrapper"); py::class_< MFIter >(m, "MFIter") .def("__repr__", [](MFIter const & mfi) { @@ -307,12 +306,10 @@ void init_MultiFab(py::module &m) { /* data access in Box index space */ .def("__iter__", [](MultiFab& mf) { - //return py::make_iterator(MFIterWrapper(mf), ValidSentinel{}); return MFIter(mf); }, // Essential: keep object alive while iterator exists py::keep_alive<0, 1>() - //py::return_value_policy::reference_internal ) ; } From f2e65152d556c59ffc041a1383f8a48b73ccc27d Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 01:00:37 -0700 Subject: [PATCH 41/43] MultiFab: test again all values make sure view & original match after modification --- tests/test_multifab.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_multifab.py b/tests/test_multifab.py index d702f36e..e3e02277 100644 --- a/tests/test_multifab.py +++ b/tests/test_multifab.py @@ -62,9 +62,16 @@ def test_mfab_loop(mfab, nghost): # all components and all indices set at once to 42 marr_np[:, :, :, :] = 42. + + # values in start & end still match? assert(marr_np[0, 0, 0, 0] == marr[bx.small_end]) assert(marr_np[-1, -1, -1, -1] == marr[bx.big_end]) + # all values for all indices match between multifab & numpy view? + for n in range(mfab.num_comp): + for i, j, k in bx: + assert(marr[i, j, k, n] == 42.) + # separate test: cupy assignment & reading # TODO From 5546509b5ae8f6a0c519cd39d653a2a607127a4d Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 01:01:59 -0700 Subject: [PATCH 42/43] Windows CI: One Build in Debug --- .github/workflows/windows.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index 5bdce2aa..10015030 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -24,7 +24,7 @@ jobs: clang: name: Clang w/o MPI runs-on: windows-latest - env: {PYAMREX_LIBDIR: build/lib/site-packages/amrex/Release/} + env: {PYAMREX_LIBDIR: build/lib/site-packages/amrex/Debug/} steps: - uses: actions/checkout@v2 - uses: seanmiddleditch/gha-setup-ninja@master @@ -38,11 +38,11 @@ jobs: -DAMReX_MPI=OFF ` -DPython_EXECUTABLE=C:\\hostedtoolcache\\windows\\python\\3.9.10\\x64\\python3.exe - cmake --build build --config Release -j 2 + cmake --build build --config Debug -j 2 python3 -m pip install -v . - name: Unit tests run: | - ctest --test-dir build -C Release --output-on-failure + ctest --test-dir build -C Debug --output-on-failure # C:\\hostedtoolcache\\windows\\python\\3.9.10\\x64\\python3.exe # C:\\hostedtoolcache\\windows\\Python\\3.10.2\\x64\\python3.exe From 04f5eb1cd2a6e34958ace56c6f8e92036b647cee Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Sat, 26 Mar 2022 01:18:48 -0700 Subject: [PATCH 43/43] Array4 __repr__: Add Type --- src/Base/Array4.cpp | 5 +++-- tests/test_array4.py | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Base/Array4.cpp b/src/Base/Array4.cpp index d460eb5e..783e9cfb 100644 --- a/src/Base/Array4.cpp +++ b/src/Base/Array4.cpp @@ -26,10 +26,11 @@ void make_Array4(py::module &m, std::string typestr) auto const array_name = std::string("Array4_").append(typestr); py::class_< Array4 >(m, array_name.c_str()) .def("__repr__", - [](Array4 const & a4) { + [typestr](Array4 const & a4) { std::stringstream s; s << a4.size(); - return ""; + return ""; } ) #if defined(AMREX_DEBUG) || defined(AMREX_BOUND_CHECK) diff --git a/tests/test_array4.py b/tests/test_array4.py index 7a94b770..31bb3008 100644 --- a/tests/test_array4.py +++ b/tests/test_array4.py @@ -23,6 +23,7 @@ def test_array4(): print(f"\nx: {x.__array_interface__} {x.dtype}") arr = amrex.Array4_double(x) print(f"arr: {arr.__array_interface__}") + print(arr) assert(arr.nComp == 1) # change original array