From 1be5b34989270fe7d43302877d4ea54360758e7a Mon Sep 17 00:00:00 2001 From: deegan Date: Mon, 9 Sep 2024 15:49:12 +0200 Subject: [PATCH] ?? --- pyphare/pyphare/core/phare_utilities.py | 12 ++- .../pharesee/hierarchy/hierarchy_utils.py | 12 ++- .../pyphare/pharesee/hierarchy/patchdata.py | 10 ++- .../samrai_hdf5_field_initializer.hpp | 61 +++++++++++---- .../initializers/samrai_hdf5_initializer.hpp | 74 ++++++++++--------- src/core/data/ndarray/ndarray_vector.hpp | 52 ++++++++----- src/core/utilities/box/box.hpp | 5 ++ .../simulator/test_samrai_restarts_parser.py | 56 ++++++-------- 8 files changed, 173 insertions(+), 109 deletions(-) diff --git a/pyphare/pyphare/core/phare_utilities.py b/pyphare/pyphare/core/phare_utilities.py index c6937cc76..29c28bb92 100644 --- a/pyphare/pyphare/core/phare_utilities.py +++ b/pyphare/pyphare/core/phare_utilities.py @@ -128,10 +128,18 @@ def is_fp32(item): return isinstance(item, float) -def assert_fp_any_all_close(a, b, atol=1e-16, rtol=0, atol_fp32=None): +def any_fp_tol(a, b, atol=1e-16, rtol=0, atol_fp32=None): if any([is_fp32(el) for el in [a, b]]): atol = atol_fp32 if atol_fp32 else atol * 1e8 - np.testing.assert_allclose(a, b, atol=atol, rtol=rtol) + return dict(atol=atol, rtol=rtol) + + +def fp_any_all_close(a, b, atol=1e-16, rtol=0, atol_fp32=None): + return np.allclose(a, b, **any_fp_tol(a, b, atol, rtol, atol_fp32)) + + +def assert_fp_any_all_close(a, b, atol=1e-16, rtol=0, atol_fp32=None): + np.testing.assert_allclose(a, b, **any_fp_tol(a, b, atol, rtol, atol_fp32)) def decode_bytes(input, errors="ignore"): diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index 2efadfc3d..8700ce4c6 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -563,6 +563,9 @@ class EqualityReport: def __bool__(self): return self.ok + def __repr__(self): + return self.reason + def hierarchy_compare(this, that): if not isinstance(this, PatchHierarchy) or not isinstance(that, PatchHierarchy): @@ -599,12 +602,7 @@ def hierarchy_compare(this, that): patch_data_cmp = patch_cmp.patch_datas[patch_data_key] if patch_data_cmp != patch_data_ref: - return EqualityReport( - False, - "data mismatch: " - + type(patch_data_cmp).__name__ - + " " - + type(patch_data_ref).__name__, - ) + msg = f"data mismatch: {patch_data_key} {type(patch_data_cmp).__name__} {type(patch_data_ref).__name__}" + return EqualityReport(False, msg) return EqualityReport(True, "OK") diff --git a/pyphare/pyphare/pharesee/hierarchy/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/patchdata.py index 717cd10ac..4b606ccdb 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patchdata.py +++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py @@ -1,6 +1,6 @@ import numpy as np -from ...core.phare_utilities import deep_copy +from ...core.phare_utilities import deep_copy, fp_any_all_close from ...core import box as boxm from ...core.box import Box @@ -80,8 +80,14 @@ def __str__(self): def __repr__(self): return self.__str__() + def compare(self, that, atol=1e-16): + return fp_any_all_close(self.dataset[:], that.dataset[:], atol) + def __eq__(self, that): - return self.field_name == that.field_name and self.dataset[:] == that.dataset[:] + return self.field_name == that.field_name and self.compare(that) + + def __ne__(self, that): + return not (self == that) def select(self, box): """ diff --git a/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp b/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp index f94d59553..369e60f21 100644 --- a/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp +++ b/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp @@ -6,15 +6,16 @@ #include #include +#include "core/def.hpp" +#include "core/logger.hpp" #include "core/data/grid/gridlayoutdefs.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" #include "core/hybrid/hybrid_quantities.hpp" #include "core/utilities/types.hpp" #include "core/data/ions/particle_initializers/particle_initializer.hpp" #include "core/data/particles/particle.hpp" #include "initializer/data_provider.hpp" #include "core/utilities/point/point.hpp" -#include "core/def.hpp" -#include "core/logger.hpp" #include "hdf5/detail/h5/h5_file.hpp" @@ -45,6 +46,7 @@ void SamraiHDF5FieldInitializer::load(Field_t& field, GridLayout const& layout) const { bool static constexpr c_ordering = false; + // using Viewer = core::NdArrayViewer; auto const local_cell = [&](auto const& box, auto const& point) { core::Point localPoint; @@ -54,47 +56,78 @@ void SamraiHDF5FieldInitializer::load(Field_t& field, return localPoint; }; - PHARE_LOG_LINE_STR("SamraiHDF5FieldInitializer::load"); + // PHARE_LOG_LINE_STR("SamraiHDF5FieldInitializer::load"); - auto const& dest_box = grow(layout.AMRBox(), GridLayout::nbrGhosts()); - auto const& overlaps = SamraiH5Interface::INSTANCE().box_intersections(dest_box); + auto const& centering = layout.centering(field.physicalQuantity()); + auto const& dest_box = layout.AMRBox(); // grow(layout.AMRBox(), GridLayout::nbrGhosts()); + auto const& overlaps = SamraiH5Interface::INSTANCE().box_intersections(dest_box); - PHARE_LOG_LINE_STR(layout.AMRBox()); + // PHARE_LOG_LINE_STR(layout.AMRBox()); for (auto const& [h5FilePtr, pdataptr] : overlaps) { auto& h5File = *h5FilePtr; auto& pdata = *pdataptr; - PHARE_LOG_LINE_STR(pdata.box); - auto const src_box = grow(pdata.box, GridLayout::nbrGhosts()); + // PHARE_LOG_LINE_STR(pdata.box); + + auto const src_box = [&]() { + auto copy = pdata.box; + // for (std::size_t i = 0; i < Field_t::dimension; ++i) + // copy.upper[i] += centering[i] == core::QtyCentering::primal ? 1 : 0; + return copy; + }(); std::vector data; std::string const fieldpath = pdata.base_path + "/" + field.name() + "##default/field_" + field.name(); h5File.file().getDataSet(fieldpath).read(data); + // PHARE_LOG_LINE_STR(data.size()); + core::Box const lcl_src_box{ core::Point{core::ConstArray()}, core::Point{ core::for_N([&](auto i) { - return static_cast(src_box.upper[i] - src_box.lower[i] + 1); + return static_cast( + src_box.upper[i] - src_box.lower[i] + (GridLayout::nbrGhosts() * 2) + + (centering[i] == core::QtyCentering::primal ? 1 : 0)); })}}; + // PHARE_LOG_LINE_STR(lcl_src_box.size()); + + assert(data.size() == lcl_src_box.size()); + auto data_view = core::make_array_view(data.data(), *lcl_src_box.shape()); auto dst_iter = dest_box.begin(); auto src_iter = src_box.begin(); for (; src_iter != src_box.end(); ++src_iter) { + // for (auto const& el : **src_iter) + // if (el < 0) + // continue; + // if (core::any(*src_iter, [](auto const& el) { el < 0; })) + // continue; // skip ghosts + if (isIn(core::Point{*src_iter}, dest_box)) { - while (*dst_iter != *src_iter) + while (*dst_iter != *src_iter and dst_iter != dest_box.end()) ++dst_iter; + + if (dst_iter == dest_box.end() or not isIn(core::Point{*dst_iter}, dest_box)) + break; + + // auto src_idx = Viewer::idx(lcl_src_box.shape(), *src_iter); + // auto dst_idx = Viewer::idx(dest_box.shape(), *dst_iter); + + // PHARE_LOG_LINE_STR(src_idx); + // PHARE_LOG_LINE_STR(dst_idx); + field(local_cell(dest_box, *dst_iter)) = data_view(local_cell(src_box, *src_iter)); - PHARE_LOG_LINE_STR(*src_iter); - PHARE_LOG_LINE_STR(*dst_iter); + // PHARE_LOG_LINE_STR(*src_iter); + // PHARE_LOG_LINE_STR(*dst_iter); - PHARE_LOG_LINE_STR(field(local_cell(dest_box, *dst_iter))); - PHARE_LOG_LINE_STR(data_view(local_cell(src_box, *src_iter))); + // PHARE_LOG_LINE_STR(field(local_cell(dest_box, *dst_iter))); + // PHARE_LOG_LINE_STR(data_view(local_cell(src_box, *src_iter))); } } } diff --git a/src/amr/data/initializers/samrai_hdf5_initializer.hpp b/src/amr/data/initializers/samrai_hdf5_initializer.hpp index a139b266a..3bcdcf655 100644 --- a/src/amr/data/initializers/samrai_hdf5_initializer.hpp +++ b/src/amr/data/initializers/samrai_hdf5_initializer.hpp @@ -87,6 +87,28 @@ class SamraiH5Interface std::unordered_map box2dataset; }; + +// struct BoxData +// { +// int dim; +// std::array lo; +// std::array hi; +// }; + +typedef struct +{ + int dim; + std::array lo; + std::array hi; +} BoxData; + +HighFive::CompoundType static box_compound_type() +{ + return {{"dim", HighFive::create_datatype()}, + {"lo", HighFive::create_datatype>()}, + {"hi", HighFive::create_datatype>()}}; +} + template struct SamraiH5Interface::SamraiHDF5File : public hdf5::h5::HighFiveFile { @@ -129,43 +151,22 @@ struct SamraiH5Interface::SamraiHDF5File : public hdf5::h5::HighFive } */ - - - struct BoxData - { - std::int32_t dim; - std::array lo; - std::array hi; - }; - - // struct BoxData - // { - // int dim; - // int lo0, lo1, lo2; - // int hi0, hi1, hi2; - // }; - - HighFive::CompoundType static box_compound_type() - { - return {{"dim", HighFive::create_datatype()}, - {"lo", HighFive::create_datatype>()}, - {"hi", HighFive::create_datatype>()}}; - } - auto getBoxFromPath(std::string const& path) const { // auto const& data = Super::read_data_set(path); - // PHARE_LOG_LINE_STR(data.size()); - std::vector boxes; - Super::file().getDataSet(path).read(boxes); - assert(boxes.size() == 1); + PHARE_LOG_LINE_STR(path); + std::vector data; + Super::file().getDataSet(path).read(data); + assert(data.size() == 1); + + // return Box_t{core::as_sized_array(data[1], data[2], data[3]), + // core::as_sized_array(data[4], data[5], data[6])}; - // auto const& bx = boxes[0]; // return Box_t{core::as_sized_array(bx.lo0, bx.lo1, bx.lo2), // core::as_sized_array(bx.hi0, bx.hi1, bx.hi2)}; - return Box_t{core::sized_array(boxes[0].lo), - core::sized_array(boxes[0].hi)}; + return Box_t{core::sized_array(data[0].lo), + core::sized_array(data[0].hi)}; } std::vector> patches; @@ -173,12 +174,17 @@ struct SamraiH5Interface::SamraiHDF5File : public hdf5::h5::HighFive - template void SamraiH5Interface::populate_from(std::string const& dir, int const& idx, int const& mpi_size) { - Box_t const mock{{0}, {99}}; + auto upstr = core::get_env("PHARE_RESTART_BUP", "199"); + int up = 199; + std::stringstream ss{upstr}; + ss >> up; + + + Box_t const mock{{0}, {up}}; for (int rank = 0; rank < mpi_size; ++rank) { auto const hdf5_filepath = getRestartFileFullPath(dir, idx, mpi_size, rank); @@ -186,7 +192,8 @@ void SamraiH5Interface::populate_from(std::string const& dir, int co for (auto const& group : h5File.scan_for_groups({"level_0000", "field_EM_B_x"})) { auto const em_path = group.substr(0, group.rfind("/")); - h5File.patches.emplace_back(mock, em_path.substr(0, em_path.rfind("/"))); + h5File.patches.emplace_back(/*h5File.getBoxFromPath(em_path + "/d_box")*/ mock, + em_path.substr(0, em_path.rfind("/"))); } } } @@ -195,5 +202,6 @@ void SamraiH5Interface::populate_from(std::string const& dir, int co } // namespace PHARE::amr +HIGHFIVE_REGISTER_TYPE(PHARE::amr::BoxData, PHARE::amr::box_compound_type) #endif diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index 7d52d1d01..72a521fe9 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -2,7 +2,9 @@ #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP #include "core/def.hpp" +#include "core/logger.hpp" #include +#include #include #include #include @@ -19,20 +21,29 @@ namespace PHARE::core template struct NdArrayViewer { + template typename Indexes, typename Index> + NO_DISCARD static std::size_t idx(NCells const& nCells, Indexes const& indexes) + + { + if constexpr (dim == 1) + return idx(nCells, indexes[0]); + + else if constexpr (dim == 2) + return idx(nCells, indexes[0], indexes[1]); + + else if constexpr (dim == 3) + return idx(nCells, indexes[0], indexes[1], indexes[2]); + } + template - NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, - Indexes const&... indexes) + NO_DISCARD static std::size_t idx(NCells const& nCells, Indexes const&... indexes) { auto params = std::forward_as_tuple(indexes...); static_assert(sizeof...(Indexes) == dim); - // static_assert((... && std::is_unsigned_v)); TODO : manage later if - // this test should be included if constexpr (dim == 1) { - auto i = std::get<0>(params); - - return data[i]; + return std::get<0>(params); } if constexpr (dim == 2) @@ -41,9 +52,9 @@ struct NdArrayViewer auto j = std::get<1>(params); if constexpr (c_ordering) - return data[j + i * nCells[1]]; + return j + i * nCells[1]; else - return data[i + j * nCells[0]]; + return i + j * nCells[0]; } if constexpr (dim == 3) @@ -53,25 +64,28 @@ struct NdArrayViewer auto k = std::get<2>(params); if constexpr (c_ordering) - return data[k + j * nCells[2] + i * nCells[1] * nCells[2]]; + return k + j * nCells[2] + i * nCells[1] * nCells[2]; else - return data[i + j * nCells[0] + k * nCells[1] * nCells[0]]; + return i + j * nCells[0] + k * nCells[1] * nCells[0]; } } + + template + NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, + Indexes const&... indexes) + { + return data[idx(nCells, indexes...)]; + } + template typename Indexes, typename Index> NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, Indexes const& indexes) { - if constexpr (dim == 1) - return at(data, nCells, indexes[0]); - - else if constexpr (dim == 2) - return at(data, nCells, indexes[0], indexes[1]); - - else if constexpr (dim == 3) - return at(data, nCells, indexes[0], indexes[1], indexes[2]); + auto const i = idx(nCells, indexes); + assert(i < product(nCells, std::size_t{1})); // debug bounds check + return data[i]; } }; diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index 0ff5d1a8a..c070e59b7 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -181,6 +181,11 @@ class box_iterator return *this; } + bool operator==(box_iterator const& other) const + { + return box_ == other.box_ and index_ == other.index_; + } + bool operator!=(box_iterator const& other) const { return box_ != other.box_ or index_ != other.index_; diff --git a/tests/simulator/test_samrai_restarts_parser.py b/tests/simulator/test_samrai_restarts_parser.py index 79bca7a3f..92ea95145 100644 --- a/tests/simulator/test_samrai_restarts_parser.py +++ b/tests/simulator/test_samrai_restarts_parser.py @@ -5,49 +5,43 @@ import datetime import unittest import numpy as np -from pathlib import Path -from datetime import timedelta -from ddt import ddt, data, unpack - -from pyphare.cpp import cpp_lib - -cpp = cpp_lib() import pyphare.pharein as ph -from pyphare.pharesee.run import Run from pyphare.simulator.simulator import Simulator from tests.simulator import SimulatorTest from tests.diagnostic import dump_all_diags from pyphare.pharesee.hierarchy.patchdata import ParticleData from pyphare.pharesee.hierarchy.fromh5 import get_all_available_quantities_from_h5 +from pyphare.pharesee.hierarchy.hierarchy_utils import hierarchy_compare -# ./build/tests/simulator/phare_outputs/restarts/test/test_restarts_1/1/1/1/00000.00400/restore.000000/nodes.0000001/proc.0000000 - -def setup_model(ppc=100): - model = ph.MaxwellianFluidModel(protons={"nbr_part_per_cell": ppc}) +def setup_model(sim, ppc=100): + model = ph.MaxwellianFluidModel( + protons={"mass": 1, "charge": 1, "nbr_part_per_cell": ppc}, + alpha={"mass": 4.0, "charge": 1, "nbr_part_per_cell": ppc}, + ) ph.ElectronModel(closure="isothermal", Te=0.12) + dump_all_diags(model.populations) return model timestep = 0.001 out = "phare_outputs/restarts/test/test_restarts_1/1/1/1/00000.00400" -simArgs = dict( +diags = "phare_outputs/restarts" + +simInitArgs = dict( + largest_patch_size=100, time_step_nbr=2, time_step=timestep, - cells=100, + cells=200, dl=0.3, init_options=dict(dir=out, mode="overwrite"), + diag_options=dict(format="phareh5", options=dict(dir=diags, mode="overwrite")), ) -def dup(dic={}): - dic.update(copy.deepcopy(simArgs)) - return dic - - def traverse_h5_for_groups_recursive(h5content: "H5Content", group, path=""): if "level_0000" in path: for key in group.attrs: @@ -70,7 +64,6 @@ def __init__(self, path): traverse_h5_for_groups_recursive(self, self.file) -@ddt class RestartsParserTest(SimulatorTest): def __init__(self, *args, **kwargs): super(RestartsParserTest, self).__init__(*args, **kwargs) @@ -85,16 +78,7 @@ def tearDown(self): def test_restart_parser(self): # h5_filepath = "phare_outputs/restarts/test/test_restarts_1/1/1/1/00000.00400/restore.000000/nodes.0000001/proc.0000000" - # h5 = H5Content(h5_filepath) - # print( - # h5.file[ - # "/PHARE_hierarchy/level_0000/level_0000-patch_0000000-block_0000000/EMAvg_B_x##default/d_ghost_box" - # ][:] - # ) - # for k in h5.data: - # print(k) - # print( # "h5.file[]", # h5.file[ @@ -102,9 +86,17 @@ def test_restart_parser(self): # ][:], # ) - sim = ph.Simulation(**dup()) - model = setup_model() - Simulator(sim).initialize() + sim = ph.Simulation(**copy.deepcopy(simInitArgs)) + model = setup_model(sim) + Simulator(sim).initialize().reset() + + datahier0 = get_all_available_quantities_from_h5(diags) + datahier1 = get_all_available_quantities_from_h5(out) + + eq = hierarchy_compare(datahier0, datahier1) + if not eq: + print(eq) + assert eq if __name__ == "__main__":