From 4ed3f65e4593b51b112d15f3b977a996e32bf288 Mon Sep 17 00:00:00 2001 From: deegan Date: Mon, 9 Mar 2026 16:00:35 +0100 Subject: [PATCH] data wrangler updates --- pyphare/pyphare/data/wrangler.py | 65 +--- pyphare/pyphare/pharesee/hierarchy/fromsim.py | 128 +++---- .../pharesee/hierarchy/hierarchy_utils.py | 27 +- pyphare/pyphare/pharesee/hierarchy/patch.py | 9 +- .../pyphare/pharesee/hierarchy/patchlevel.py | 18 +- src/amr/physical_models/mhd_model.hpp | 7 +- src/core/utilities/mpi_utils.cpp | 2 +- src/core/utilities/mpi_utils.hpp | 35 +- src/core/utilities/span.hpp | 37 +- src/python3/cpp_etc.cpp | 42 ++- src/python3/cpp_simulator.hpp | 48 ++- src/python3/data_wrangler.hpp | 104 +++--- src/python3/patch_data.hpp | 48 +-- src/python3/patch_level.hpp | 331 +++++------------- src/python3/pybind_def.hpp | 58 ++- tests/simulator/CMakeLists.txt | 3 +- tests/simulator/refined_particle_nbr.py | 132 ------- tests/simulator/test_data_wrangler.py | 57 +++ 18 files changed, 496 insertions(+), 655 deletions(-) delete mode 100644 tests/simulator/refined_particle_nbr.py create mode 100644 tests/simulator/test_data_wrangler.py diff --git a/pyphare/pyphare/data/wrangler.py b/pyphare/pyphare/data/wrangler.py index cf362e186..54ffac3c7 100644 --- a/pyphare/pyphare/data/wrangler.py +++ b/pyphare/pyphare/data/wrangler.py @@ -12,6 +12,8 @@ def __init__(self, simulator): self.cpp = getattr(cpp.cpp_lib(sim), "DataWrangler")( simulator.cpp_sim, simulator.cpp_hier ) + # todo mhd + self.modelsPerLevel = ["Hybrid" for lvl in range(self.cpp.getNumberOfLevels())] def kill(self): del self.cpp @@ -19,55 +21,18 @@ def kill(self): def getNumberOfLevels(self): return self.cpp.getNumberOfLevels() - def getPatchLevel(self, lvl): - return self.cpp.getPatchLevel(lvl) - - def _lvl0FullContiguous(self, input, is_primal=True): - return self.cpp.sync_merge(input, is_primal) - - def lvl0IonDensity(self): - return self._lvl0FullContiguous(self.getPatchLevel(0).getDensity()) - - def lvl0BulkVelocity(self): - return { - xyz: self._lvl0FullContiguous(bv) - for xyz, bv in self.getPatchLevel(0).getBulkVelocity().items() - } - - def lvl0PopDensity(self): - return { - pop: self._lvl0FullContiguous(density) - for pop, density in self.getPatchLevel(0).getPopDensities().items() - } + def getMHDPatchLevel(self, lvl): + return self.cpp.getMHDPatchLevel(lvl) - def lvl0PopFluxes(self): - return { - pop: {xyz: self._lvl0FullContiguous(data) for xyz, data in flux.items()} - for pop, flux in self.getPatchLevel(0).getPopFluxes().items() - } + def getHybridPatchLevel(self, lvl): + return self.cpp.getHybridPatchLevel(lvl) - def extract_is_primal_key_from(self, em_xyz): - """extract "Ex" from "EM_E_x" """ - return "".join(em_xyz.split("_"))[2:] - - def lvl0EM(self): - return { - em: { - em_xyz: self.cpp.sync_merge( - data, - gridlayout.yee_element_is_primal( - self.extract_is_primal_key_from(em_xyz) - ), - ) - for em_xyz, data in xyz_map.items() - } - for em, xyz_map in self.getPatchLevel(0).getEM().items() - } - - -# for pop, particles in dw.getPatchLevel(0).getParticles().items(): -# print("pop :", pop) -# for key, patches in particles.items(): -# print("\tkey :", key) -# for patch in patches: -# print("\t\t", patch.patchID, "size:", patch.data.size()) + def getPatchLevel(self, lvl): + lvlModel = self.modelsPerLevel[lvl] + assert lvlModel in ["Hybrid", "MHD"] + if lvlModel == "MHD": + return self.getMHDPatchLevel(lvl) + return self.getHybridPatchLevel(lvl) + + def sync(self, patch_datas): + return self.cpp.sync(patch_datas) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromsim.py b/pyphare/pyphare/pharesee/hierarchy/fromsim.py index aa13ea2cc..fe69a5bc2 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromsim.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromsim.py @@ -1,3 +1,4 @@ +from pyphare import cpp from .hierarchy_utils import isFieldQty, field_qties, quantidic, refinement_ratio from .patchdata import FieldData, ParticleData from .patch import Patch @@ -10,14 +11,25 @@ import numpy as np -def hierarchy_from_sim(simulator, qty, pop=""): +def patch_gridlayout(patch, dl, interp_order): + return GridLayout( + Box(patch.lower, patch.upper), patch.origin, dl, interp_order=interp_order + ) + + +def hierarchy_from_sim(simulator, qty, pop="", hier=None, sync=False): + """ + sync==True will copy all data to rank 0! + leaving other ranks with an empty hierarch! + """ + dw = simulator.data_wrangler() nbr_levels = dw.getNumberOfLevels() patch_levels = {} - root_cell_width = simulator.cell_width() + root_cell_width = np.asarray(simulator.cell_width()) domain_box = Box([0] * len(root_cell_width), simulator.domain_box()) - assert len(domain_box.ndim) == len(simulator.domain_box().ndim) + assert domain_box.ndim == len(simulator.domain_box()) for ilvl in range(nbr_levels): lvl_cell_width = root_cell_width / refinement_ratio**ilvl @@ -26,98 +38,46 @@ def hierarchy_from_sim(simulator, qty, pop=""): getters = quantidic(ilvl, dw) if isFieldQty(qty): - wpatches = getters[qty]() - for patch in wpatches: - patch_datas = {} - lower = patch.lower - upper = patch.upper - origin = patch.origin - layout = GridLayout( - Box(lower, upper), - origin, - lvl_cell_width, - interp_order=simulator.interporder(), + patch_datas = getters[qty]() + + if sync: + patch_datas = dw.sync(patch_datas) + if cpp.mpi_rank() > 0: + continue + + for patch_data in patch_datas: + layout = patch_gridlayout( + patch_data, lvl_cell_width, simulator.interp_order() + ) + patches[ilvl].append( + Patch({qty: FieldData(layout, field_qties[qty], patch_data.data)}) ) - pdata = FieldData(layout, field_qties[qty], patch.data) - patch_datas[qty] = pdata - patches[ilvl].append(Patch(patch_datas)) - elif qty == "particles": + elif qty == "particles": # domain only! + if sync: + raise ValueError("sync not supported for particles") if pop == "": raise ValueError("must specify pop argument for particles") - # here the getter returns a dict like this - # {'protons': {'patchGhost': [, - # ], - # 'domain': [, - # ]}} - - # domain particles are assumed to always be here - # but patchGhost and levelGhost may not be, depending on the level - - populationdict = getters[qty](pop)[pop] - - dom_dw_patches = populationdict["domain"] - for patch in dom_dw_patches: - patch_datas = {} - - lower = patch.lower - upper = patch.upper - origin = patch.origin - layout = GridLayout( - Box(lower, upper), - origin, - lvl_cell_width, - interp_order=simulator.interp_order(), - ) - v = np.asarray(patch.data.v).reshape(int(len(patch.data.v) / 3), 3) - - domain_particles = Particles( - icells=np.asarray(patch.data.iCell), - deltas=np.asarray(patch.data.delta), - v=v, - weights=np.asarray(patch.data.weight), - charges=np.asarray(patch.data.charge), + + for patch_data in getters[qty](pop): + layout = patch_gridlayout( + patch_data, lvl_cell_width, simulator.interp_order() ) - patch_datas[pop + "_particles"] = ParticleData( - layout, domain_particles, pop + patches[ilvl].append( # ParticleData is SoA COPY! + Patch({qty: FieldData(layout, "tags", patch_data.data)}) ) - patches[ilvl].append(Patch(patch_datas)) - - # ok now let's add the patchGhost if present - # note that patchGhost patches may not be the same list as the - # domain patches... since not all patches may not have patchGhost while they do have - # domain... while looping on the patchGhost items, we need to search in - # the already created patches which one to which add the patchGhost particles - - for ghostParticles in ["levelGhost"]: - if ghostParticles in populationdict: - for dwpatch in populationdict[ghostParticles]: - v = np.asarray(dwpatch.data.v) - s = v.size - v = v[:].reshape(int(s / 3), 3) - - patchGhost_part = Particles( - icells=np.asarray(dwpatch.data.iCell), - deltas=np.asarray(dwpatch.data.delta), - v=v, - weights=np.asarray(dwpatch.data.weight), - charges=np.asarray(dwpatch.data.charge), - ) - - box = Box(dwpatch.lower, dwpatch.upper) - - # now search which of the already created patches has the same box - # once found we add the new particles to the ones already present - - patch = [p for p in patches[ilvl] if p.box == box][0] - patch.patch_datas[pop + "_particles"].dataset.add( - patchGhost_part - ) else: raise ValueError("{} is not a valid quantity".format(qty)) patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) + if hier: + for lvl_nbr, level in hier.levels(hier.times()[0]).items(): + new_level = patch_levels[lvl_nbr] + for ip, patch in enumerate(level.patches): + patch.patch_datas = {**patch.patch_datas, **new_level[ip].patch_datas} + + return hier return PatchHierarchy(patch_levels, domain_box, time=simulator.currentTime()) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index d3856873d..5d8ecd5c6 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -204,19 +204,20 @@ def quantidic(ilvl, wrangler): pl = wrangler.getPatchLevel(ilvl) return { - "density": pl.getDensity, - "bulkVelocity_x": pl.getVix, - "bulkVelocity_y": pl.getViy, - "bulkVelocity_z": pl.getViz, - "EM_B_x": pl.getBx, - "EM_B_y": pl.getBy, - "EM_B_z": pl.getBz, - "EM_E_x": pl.getEx, - "EM_E_y": pl.getEy, - "EM_E_z": pl.getEz, - "flux_x": pl.getFx, - "flux_y": pl.getFy, - "flux_z": pl.getFz, + "density": pl.getNi, + "Vi": pl.getVi, + "EM_B_x": lambda: pl.getB("x"), + "EM_B_y": lambda: pl.getB("y"), + "EM_B_z": lambda: pl.getB("z"), + "EM_E_x": lambda: pl.getE("x"), + "EM_E_y": lambda: pl.getE("y"), + "EM_E_z": lambda: pl.getE("z"), + "bulkVelocity_x": lambda: pl.getVi("x"), + "bulkVelocity_y": lambda: pl.getVi("y"), + "bulkVelocity_z": lambda: pl.getVi("z"), + "flux_x": lambda pop: pl.getFlux("x", pop), + "flux_y": lambda pop: pl.getFlux("y", pop), + "flux_z": lambda pop: pl.getFlux("z", pop), "particles": pl.getParticles, } diff --git a/pyphare/pyphare/pharesee/hierarchy/patch.py b/pyphare/pyphare/pharesee/hierarchy/patch.py index 0ac18f3d8..488af1bf3 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patch.py +++ b/pyphare/pyphare/pharesee/hierarchy/patch.py @@ -37,8 +37,15 @@ def __str__(self): def __repr__(self): return self.__str__() + def __iter__(self): + return self.patch_datas.items().__iter__() + def __getitem__(self, key): - return self.patch_datas[key] + try: + return self.patch_datas[key] + except KeyError as e: + print(key, "not in", self.patch_datas.keys()) + raise e def copy(self): """does not copy patchdatas.datasets (see class PatchData)""" diff --git a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py index 72fee0d5d..00e94753a 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py +++ b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py @@ -1,3 +1,8 @@ +# +# +# + + class PatchLevel: """is a collection of patches""" @@ -5,11 +10,18 @@ def __init__(self, lvl_nbr, patches): self.level_number = lvl_nbr self.patches = patches - def __iter__(self): - return self.patches.__iter__() - def level_range(self): name = list(self.patches[0].patch_datas.keys())[0] return min([patch.patch_datas[name].x.min() for patch in self.patches]), max( [patch.patch_datas[name].x.max() for patch in self.patches] ) + + def __getitem__(self, idx): + return self.patches[idx] + + def __iter__(self): + return self.patches.__iter__() + + @property + def cell_width(self): + return self.patches[0].layout.dl diff --git a/src/amr/physical_models/mhd_model.hpp b/src/amr/physical_models/mhd_model.hpp index 1f155f7ee..a8d0133c0 100644 --- a/src/amr/physical_models/mhd_model.hpp +++ b/src/amr/physical_models/mhd_model.hpp @@ -23,9 +23,10 @@ namespace solver class MHDModel : public IPhysicalModel { public: - using patch_t = typename AMR_Types::patch_t; - using level_t = typename AMR_Types::level_t; - using Interface = IPhysicalModel; + using gridlayout_type = GridLayoutT; + using patch_t = AMR_Types::patch_t; + using level_t = AMR_Types::level_t; + using Interface = IPhysicalModel; static constexpr std::string_view model_type_name = "MHDModel"; static inline std::string const model_name{model_type_name}; diff --git a/src/core/utilities/mpi_utils.cpp b/src/core/utilities/mpi_utils.cpp index dae6f2b49..7779b0252 100644 --- a/src/core/utilities/mpi_utils.cpp +++ b/src/core/utilities/mpi_utils.cpp @@ -42,7 +42,7 @@ void barrier() -std::string date_time(std::string format) +std::string date_time(std::string const& format) { std::time_t t = std::time(NULL); char buffer[80]; diff --git a/src/core/utilities/mpi_utils.hpp b/src/core/utilities/mpi_utils.hpp index c093e34d3..01334b574 100644 --- a/src/core/utilities/mpi_utils.hpp +++ b/src/core/utilities/mpi_utils.hpp @@ -26,7 +26,7 @@ NO_DISCARD int rank(); void barrier(); -NO_DISCARD std::string date_time(std::string format = "%Y-%m-%d-%H:%M:%S"); +NO_DISCARD std::string date_time(std::string const& format = "%Y-%m-%d-%H:%M:%S"); NO_DISCARD std::int64_t unix_timestamp_now(); @@ -60,6 +60,11 @@ NO_DISCARD auto mpi_type_for() // don't return anything = compile failure if tried to use this function } +auto mpi_type_for(auto const& val) +{ + return mpi_type_for>(); +} + template auto all_get_from(int const& rank_, Fn&& fn, Args&&... args) @@ -142,10 +147,11 @@ void _collect_vector(SendBuff const& sendBuff, RcvBuff& rcvBuff, std::vector -NO_DISCARD std::vector collectVector(Vector const& sendBuff, int mpi_size = 0) + +template> +NO_DISCARD Return_t collectVector(Vector const& sendBuff, int mpi_size = 0) { - using Data = typename Vector::value_type; + using Data = Vector::value_type; if (mpi_size == 0) MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); @@ -156,10 +162,10 @@ NO_DISCARD std::vector collectVector(Vector const& sendBuff, int mpi_siz _collect_vector(sendBuff, rcvBuff, perMPISize, displs, mpi_size); std::size_t offset = 0; - std::vector collected; - for (int i = 0; i < mpi_size; i++) + Return_t collected; + for (int i = 0; i < mpi_size; ++i) { - if (perMPISize[i] == 0) + if (perMPISize[i] == 0) // NO OUT OF BOUNDS ON LAST RANK collected.emplace_back(); else collected.emplace_back(&rcvBuff[offset], &rcvBuff[offset] + perMPISize[i]); @@ -191,13 +197,12 @@ NO_DISCARD auto collectArrays(std::array const& arr, int mpi_size) if (mpi_size == 0) MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); - std::size_t maxMPISize = arr.size(); - std::vector datas(maxMPISize * mpi_size); - _collect(arr.data(), datas, arr.size(), maxMPISize); + std::vector datas(size * mpi_size); + _collect(arr.data(), datas, size, size); std::vector values(mpi_size); for (int i = 0; i < mpi_size; i++) - std::memcpy(&values[i], &datas[maxMPISize * i], maxMPISize); + std::memcpy(&values[i], &datas[size * i], size * sizeof(T)); return values; } @@ -213,6 +218,14 @@ NO_DISCARD SpanSet collect_raw(Vector const& d } +template +NO_DISCARD auto collect(Span const& sendBuff, int mpi_size = 0) +{ + using V = Span::value_type; + return collectVector, std::vector>>(sendBuff, mpi_size); +} + + template NO_DISCARD std::vector collect(Data const& data, int mpi_size) { diff --git a/src/core/utilities/span.hpp b/src/core/utilities/span.hpp index dab49c849..6e1b80cf7 100644 --- a/src/core/utilities/span.hpp +++ b/src/core/utilities/span.hpp @@ -3,6 +3,7 @@ #ifndef PHARE_CORE_UTILITIES_SPAN_HPP #define PHARE_CORE_UTILITIES_SPAN_HPP + #include "core/def.hpp" #include "core/utilities/types.hpp" @@ -21,20 +22,35 @@ concept Spannable = requires(T t) { template - struct Span { - using value_type = T; + using value_type = std::decay_t; + + Span(T* ptr_ = nullptr, SIZE const s_ = 0) + : ptr{ptr_} + , s{s_} + { + } + + Span(Span&&) = default; + Span(Span const&) = default; + Span& operator=(Span&&) = default; + Span& operator=(Span const&) = default; NO_DISCARD auto& operator[](SIZE i) { return ptr[i]; } NO_DISCARD auto& operator[](SIZE i) const { return ptr[i]; } - NO_DISCARD T const* const& data() const { return ptr; } - NO_DISCARD T const* const& begin() const { return ptr; } - NO_DISCARD T* end() const { return ptr + s; } + // NO_DISCARD T const* cdata() const { return ptr; } + NO_DISCARD auto data() const { return ptr; } + NO_DISCARD auto data() { return ptr; } + NO_DISCARD auto begin() { return ptr; } + NO_DISCARD auto begin() const { return ptr; } + NO_DISCARD auto end() { return ptr + s; } + NO_DISCARD auto end() const { return ptr + s; } NO_DISCARD SIZE const& size() const { return s; } + NO_DISCARD auto size_address() { return &s; } - T const* ptr = nullptr; - SIZE s = 0; + T* ptr = nullptr; + SIZE s = 0; }; @@ -88,7 +104,12 @@ struct SpanSet { } - NO_DISCARD Span operator[](SIZE i) const + + NO_DISCARD Span operator[](SIZE i) + { + return {this->vec.data() + displs[i], this->sizes[i]}; + } + NO_DISCARD Span operator[](SIZE i) const { return {this->vec.data() + displs[i], this->sizes[i]}; } diff --git a/src/python3/cpp_etc.cpp b/src/python3/cpp_etc.cpp index 6c13439aa..e04834ea0 100644 --- a/src/python3/cpp_etc.cpp +++ b/src/python3/cpp_etc.cpp @@ -17,11 +17,33 @@ #include "hdf5/detail/h5/h5_file.hpp" #endif + +#include + namespace py = pybind11; namespace PHARE::pydata { +template +void declareParticles(py::module& m) +{ + using Particle = core::Particle; + std::string name = "Particle_" + std::to_string(dim) + "D"; + py::class_>(m, name.c_str()) + .def_readonly("iCell", &Particle::iCell) + .def_readonly("delta", &Particle::delta) + .def_readonly("v", &Particle::v); + + using ParticleArray = core::ParticleArray; + name = "ParticleArray_" + std::to_string(dim) + "D"; + py::class_>(m, name.c_str()) + .def("size", &ParticleArray::size) + .def("__getitem__", + [](ParticleArray& self, std::size_t const idx) -> auto& { return self[idx]; }); +} + + template void declarePatchData(py::module& m, std::string key) { @@ -88,7 +110,7 @@ PYBIND11_MODULE(cpp_etc, m) .def(py::init<>()) .def("reset", &SamraiLifeCycle::reset); - py::class_>(m, "AMRHierarchy"); + py::class_(m, "AMRHierarchy"); m.def("make_hierarchy", []() { return PHARE::amr::Hierarchy::make(); }); m.def("makePyArrayWrapper", makePyArrayWrapper); @@ -134,9 +156,21 @@ PYBIND11_MODULE(cpp_etc, m) declareDim<2>(m); declareDim<3>(m); - declarePatchData, 1>(m, "PatchDataVectorDouble_1D"); - declarePatchData, 2>(m, "PatchDataVectorDouble_2D"); - declarePatchData, 3>(m, "PatchDataVectorDouble_3D"); + declarePatchData, 1>(m, "PatchPyArrayDouble_1D"); + declarePatchData, 2>(m, "PatchPyArrayDouble_2D"); + declarePatchData, 3>(m, "PatchPyArrayDouble_3D"); + + declareParticles<1>(m); + declareParticles<2>(m); + declareParticles<3>(m); + + declarePatchData, 1>(m, "PatchDataParticleArray_1D"); + declarePatchData, 2>(m, "PatchDataParticleArray_2D"); + declarePatchData, 3>(m, "PatchDataParticleArray_3D"); + + declarePatchData*, 1>(m, "PatchDataParticleArrayPtr_1D"); + declarePatchData*, 2>(m, "PatchDataParticleArrayPtr_2D"); + declarePatchData*, 3>(m, "PatchDataParticleArrayPtr_3D"); } } // namespace PHARE::pydata diff --git a/src/python3/cpp_simulator.hpp b/src/python3/cpp_simulator.hpp index fff92aa11..bc38cbfd7 100644 --- a/src/python3/cpp_simulator.hpp +++ b/src/python3/cpp_simulator.hpp @@ -59,33 +59,31 @@ void inline declare_etc(py::module& m) py::class_(m, name.c_str()) .def(py::init const&, std::shared_ptr const&>()) .def(py::init const&, std::shared_ptr const&>()) - .def("sync_merge", &DW::sync_merge) - .def("getPatchLevel", &DW::getPatchLevel) + .def("sync", &DW::sync) + .def("getMHDPatchLevel", &DW::getMHDPatchLevel) + .def("getHybridPatchLevel", &DW::getHybridPatchLevel) .def("getNumberOfLevels", &DW::getNumberOfLevels); - using PL = PatchLevel; - name = "PatchLevel"; - py::class_(m, name.c_str()) - .def("getEM", &PL::getEM) - .def("getE", &PL::getE) - .def("getB", &PL::getB) - .def("getBx", &PL::getBx) - .def("getBy", &PL::getBy) - .def("getBz", &PL::getBz) - .def("getEx", &PL::getEx) - .def("getEy", &PL::getEy) - .def("getEz", &PL::getEz) - .def("getVix", &PL::getVix) - .def("getViy", &PL::getViy) - .def("getViz", &PL::getViz) - .def("getDensity", &PL::getDensity) - .def("getBulkVelocity", &PL::getBulkVelocity) - .def("getPopDensities", &PL::getPopDensities) - .def("getPopFluxes", &PL::getPopFlux) - .def("getFx", &PL::getFx) - .def("getFy", &PL::getFy) - .def("getFz", &PL::getFz) - .def("getParticles", &PL::getParticles, py::arg("userPopName") = "all"); + + using HybPL = PatchLevel; + name = "HybridPatchLevel"; + if constexpr (core::defaultNbrRefinedParts(opts.dimension, opts.interp_order) + == opts.nbRefinedPart) // register once! + py::class_(m, name.c_str()) + .def("getB", &HybPL::getB, py::arg("component")) + .def("getE", &HybPL::getE, py::arg("component")) + .def("getVi", &HybPL::getVi, py::arg("component")) + .def("getN", &HybPL::getN, py::arg("pop_name")) + .def("getNi", &HybPL::getNi) + .def("getFlux", &HybPL::getFlux, py::arg("component"), py::arg("pop_name")) + .def("getParticles", &HybPL::getParticles, py::arg("pop_name")); + + + using MHDPL = PatchLevel; + name = "MHDPatchLevel"; + if constexpr (core::defaultNbrRefinedParts(opts.dimension, opts.interp_order) + == opts.nbRefinedPart) + py::class_(m, name.c_str()); using _Splitter = PHARE::amr::Splitter, core::InterpConst, diff --git a/src/python3/data_wrangler.hpp b/src/python3/data_wrangler.hpp index 05628bdc7..6de230a19 100644 --- a/src/python3/data_wrangler.hpp +++ b/src/python3/data_wrangler.hpp @@ -1,28 +1,34 @@ #ifndef PHARE_PYTHON_DATA_WRANGLER_HPP #define PHARE_PYTHON_DATA_WRANGLER_HPP + #include "core/utilities/mpi_utils.hpp" #include "core/utilities/point/point.hpp" #include "core/utilities/meta/meta_utilities.hpp" #include "amr/wrappers/hierarchy.hpp" +#include "core/utilities/types.hpp" #include "initializer/data_provider.hpp" #include "python3/patch_data.hpp" #include "python3/patch_level.hpp" +#include "python3/pybind_def.hpp" #include "simulator/simulator.hpp" #include "dict.hpp" + #include #include #include #include #include -#include #include +#include + + namespace PHARE::pydata { @@ -71,6 +77,8 @@ class __attribute__((visibility("hidden"))) DataWrangler using Simulator = PHARE::Simulator; using HybridModel = Simulator::HybridModel; + using MHDModel = Simulator::MHDModel; + DataWrangler(std::shared_ptr const& simulator, @@ -90,91 +98,63 @@ class __attribute__((visibility("hidden"))) DataWrangler auto getNumberOfLevels() const { return hierarchy_->getNumberOfLevels(); } - auto getPatchLevel(size_t lvl) + auto getMHDPatchLevel(size_t lvl) { - return PatchLevel{*hierarchy_, *simulator_.getHybridModel(), lvl}; + return PatchLevel{*hierarchy_, *simulator_.getMHDModel(), lvl}; } - auto sort_merge_1d(std::vector, dimension>> const&& input, - bool shared_patch_border = false) + auto getHybridPatchLevel(size_t lvl) { - std::vector, dimension> const*>> sorted; - for (auto const& data : input) - sorted.emplace_back(core::Point::fromString(data.origin)[0], &data); - std::sort(sorted.begin(), sorted.end(), [](auto& a, auto& b) { return a.first < b.first; }); - std::vector ret; - for (size_t i = 0; i < sorted.size(); i++) - { // skip empty patches in case of unequal patches across MPI domains - if (!sorted[i].second->data.size()) - continue; - auto& data = sorted[i].second->data; - auto& ghosts = sorted[i].second->nGhosts; - auto end = ghosts; - // primal nodes share a cell wall when patches touch so drop duplicate value if so - if (shared_patch_border) - end = i == sorted.size() - 1 ? end : end + 1; - ret.insert(std::end(ret), std::begin(data) + ghosts, std::end(data) - end); - } - return ret; + return PatchLevel{*hierarchy_, *simulator_.getHybridModel(), lvl}; } - auto sync(std::vector, dimension>> const& input) + + auto sync(std::vector, dimension>> const& input) { - int mpi_size = core::mpi::size(); - std::vector, dimension>> collected; + // collect all data on rank 0! - auto reinterpret_array = [&](auto& py_array) { - return reinterpret_cast&>( - *static_cast(py_array.request().ptr)); - }; + auto const mpi_size = core::mpi::size(); + std::vector, dimension>> collected; - auto collect = [&](PatchData, dimension> const& patch_data) { + auto const collect = [&](PatchData, dimension> const& patch_data) { auto patchIDs = core::mpi::collect(patch_data.patchID, mpi_size); - auto origins = core::mpi::collect(patch_data.origin, mpi_size); - auto lower = core::mpi::collect_raw(makeSpan(patch_data.lower), mpi_size); - auto upper = core::mpi::collect_raw(makeSpan(patch_data.upper), mpi_size); + auto shapes = core::mpi::collectArrays(shape(patch_data.data), mpi_size); + auto origins = core::mpi::collect(makeSpan(patch_data.origin), mpi_size); + auto lower = core::mpi::collect(makeSpan(patch_data.lower), mpi_size); + auto upper = core::mpi::collect(makeSpan(patch_data.upper), mpi_size); auto ghosts = core::mpi::collect(patch_data.nGhosts, mpi_size); - auto datas = core::mpi::collect(patch_data.data, mpi_size); - - for (int i = 0; i < mpi_size; i++) - { - auto& data = collected.emplace_back(); - setPatchData(data, patchIDs[i], origins[i], lower[i], upper[i]); - data.nGhosts = ghosts[i]; - data.data = std::move(datas[i]); - } + auto datas = core::mpi::collect(makeSpan(patch_data.data), mpi_size); + + if (core::mpi::rank() == 0) + for (int i = 0; i < mpi_size; ++i) + { + if (datas[i].size() == 0) // missing patch on rank + continue; + + auto& data = collected.emplace_back(shapes[i], strides_from(shapes[i])); + auto const span = makeSpan(data.data); + std::memcpy(span.data(), datas[i].data(), span.size() * sizeof(double)); + setPatchData(data, patchIDs[i], origins[i], lower[i], upper[i]); + data.nGhosts = ghosts[i]; + } }; - std::size_t max = core::mpi::max(input.size(), mpi_size); + auto const max = core::mpi::max(input.size()); - PatchData, dimension> empty; + PatchData, dimension> empty{core::ConstArray()}; + + for (std::size_t i = 0; i < max; ++i) + collect(i < input.size() ? input[i] : empty); - for (size_t i = 0; i < max; i++) - { - if (i < input.size()) - collect(input[i]); - else - collect(empty); - } return collected; } - auto sync_merge(std::vector, dimension>> const& input, - [[maybe_unused]] bool primal) - { - if constexpr (dimension == 1) - return sort_merge_1d(sync(input), primal); - - throw std::runtime_error("Not handled for >1 dim"); - } private: Simulator& simulator_; std::shared_ptr hierarchy_; - - static Simulator& cast_simulator(std::shared_ptr const& simulator) { using SimulatorCaster = SimulatorCaster; diff --git a/src/python3/patch_data.hpp b/src/python3/patch_data.hpp index e87477e29..440587f39 100644 --- a/src/python3/patch_data.hpp +++ b/src/python3/patch_data.hpp @@ -1,8 +1,12 @@ #ifndef PHARE_PYTHON_PATCH_DATA_HPP #define PHARE_PYTHON_PATCH_DATA_HPP + #include "pybind_def.hpp" +#include "core/utilities/point/point.hpp" + +#include #include #include #include @@ -10,16 +14,12 @@ namespace PHARE::pydata { +namespace py = pybind11; + template struct __attribute__((visibility("hidden"))) PatchData { static auto constexpr dimension = dim; - std::string patchID; - std::string origin; - py_array_t lower{dim}; - py_array_t upper{dim}; - std::size_t nGhosts; - Data data; PatchData() = default; @@ -29,36 +29,44 @@ struct __attribute__((visibility("hidden"))) PatchData : data{std::forward(args)...} { } + + Data data; + std::string patchID; + py_array_t origin{dim}; + py_array_t lower{dim}; + py_array_t upper{dim}; + std::size_t nGhosts; }; -template -void setPatchData(PatchData& data, std::string const patchID, std::string const origin, - Container lower, Container upper) +template +void setPatchData(PatchData& data, std::string const patchID, auto const origin, auto const lower, + auto const upper) { - constexpr std::size_t bytes = PatchData::dimension * sizeof(size_t); - std::memcpy(data.lower.request().ptr, lower.data(), bytes); - std::memcpy(data.upper.request().ptr, upper.data(), bytes); + std::memcpy(data.lower.request().ptr, lower.data(), PatchData::dimension * sizeof(int)); + std::memcpy(data.upper.request().ptr, upper.data(), PatchData::dimension * sizeof(int)); + std::memcpy(data.origin.request().ptr, origin.data(), PatchData::dimension * sizeof(double)); data.patchID = patchID; - data.origin = origin; } template -void setPatchDataFromGrid(PatchData& pdata, GridLayout& grid, std::string patchID) +void setPatchDataFromGrid(PatchData& pdata, GridLayout& grid, std::string const& patchID) { - setPatchData(pdata, patchID, grid.origin().str(), - grid.AMRBox().lower.template toArray(), - grid.AMRBox().upper.template toArray()); + setPatchData(pdata, patchID, *grid.origin(), *grid.AMRBox().lower, *grid.AMRBox().upper); } + template -void setPatchDataFromField(PatchData& pdata, Field const& field, GridLayout& grid, - std::string patchID) +void setPyPatchDataFromField(PatchData& pdata, Field const& field, GridLayout& grid, + std::string patchID) { + auto constexpr dimension = PatchData::dimension; + static_assert(dimension >= 1 and dimension <= 3); + setPatchDataFromGrid(pdata, grid, patchID); pdata.nGhosts = static_cast( GridLayout::nbrGhosts(GridLayout::centering(field.physicalQuantity())[0])); - pdata.data.assign(field.data(), field.data() + field.size()); + pdata.data = field_as_memory_view(field); } diff --git a/src/python3/patch_level.hpp b/src/python3/patch_level.hpp index 44d141f76..61021ac5b 100644 --- a/src/python3/patch_level.hpp +++ b/src/python3/patch_level.hpp @@ -1,306 +1,169 @@ #ifndef PHARE_PYTHON_PATCH_LEVEL_HPP #define PHARE_PYTHON_PATCH_LEVEL_HPP -#include "phare_solver.hpp" + + +#include "amr/wrappers/hierarchy.hpp" +#include "amr/physical_models/mhd_model.hpp" +#include "amr/physical_models/hybrid_model.hpp" + #include "python3/patch_data.hpp" + + #include #include #include + namespace PHARE::pydata { -template -class __attribute__((visibility("hidden"))) PatchLevel + + +template +class PatchLevel; + +template +class AnyPatchLevel { -public: - static constexpr std::size_t dimension = opts.dimension; + using GridLayout = Model_t::gridlayout_type; - using PHARESolverTypes = solver::PHARE_Types; - using HybridModel = PHARESolverTypes::HybridModel_t; - using GridLayout = HybridModel::gridlayout_type; + static constexpr std::size_t dimension = GridLayout::dimension; - PatchLevel(amr::Hierarchy& hierarchy, HybridModel& model, std::size_t lvl) - : lvl_(lvl) - , hierarchy_{hierarchy} - , model_{model} +public: + AnyPatchLevel(amr::Hierarchy& hierarchy, Model_t& model, std::size_t lvl) + : lvl(lvl) + , hierarchy{hierarchy} + , model{model} { } - auto getDensity() - { - std::vector, dimension>> patchDatas; - auto& ions = model_.state.ions; - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - setPatchDataFromField(patchDatas.emplace_back(), ions.chargeDensity(), grid, patchID); +protected: + auto getField(auto& field) + { + std::vector, dimension>> patchDatas; + auto visit = [&](auto& grid, auto patchID, auto /*ilvl*/) { + setPyPatchDataFromField(patchDatas.emplace_back(), field, grid, patchID); }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); - + amr::visitLevel(*hierarchy.getPatchLevel(lvl), rm, visit, field); return patchDatas; } - - auto getPopDensities() + auto getVecFieldComponent(auto& vecfield, auto const component) { - using Inner = decltype(getDensity()); - - std::unordered_map pop_data; - auto& ions = model_.state.ions; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - for (auto const& pop : ions) - { - if (!pop_data.count(pop.name())) - pop_data.emplace(pop.name(), Inner()); - - setPatchDataFromField(pop_data.at(pop.name()).emplace_back(), pop.chargeDensity(), - grid, patchID); - } + std::vector, dimension>> patchDatas; + auto visit = [&](auto& grid, auto patchID, auto /*ilvl*/) { + setPyPatchDataFromField(patchDatas.emplace_back(), vecfield(component), grid, patchID); }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); - - return pop_data; + amr::visitLevel(*hierarchy.getPatchLevel(lvl), rm, visit, vecfield); + return patchDatas; } - template - auto getVecFields(VecField& vecField, Map& container, GridLayout& grid, std::string patchID, - std::string outer_key) + auto static vecfield_component(std::string const& component) { - if (!container.count(outer_key)) - container.emplace(outer_key, typename Map::mapped_type()); - - auto& inner = container.at(outer_key); - - for (auto& [id, type] : core::Components::componentMap()) - { - auto& field = vecField.getComponent(type); - - if (!inner.count(field.name())) - inner.emplace(field.name(), decltype(getDensity())()); - - setPatchDataFromField(inner.at(field.name()).emplace_back(), field, grid, patchID); - } + return core::Components::componentMap().at(component); } - auto getBulkVelocity() - { - decltype(getPopDensities()) bulkV; - - auto& ions = model_.state.ions; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - for (auto& [id, type] : core::Components::componentMap()) - { - auto& field = ions.velocity().getComponent(type); +protected: + std::size_t lvl; + amr::Hierarchy& hierarchy; + Model_t& model; + Model_t::resources_manager_type& rm = *model.resourcesManager; +}; - if (!bulkV.count(field.name())) - bulkV.emplace(field.name(), decltype(getDensity())()); - setPatchDataFromField(bulkV.at(field.name()).emplace_back(), field, grid, patchID); - } - }; +template +class __attribute__((visibility("hidden"))) +PatchLevel>> : AnyPatchLevel +{ + using Super = AnyPatchLevel; - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); +public: + using Model_t = Model; + using GridLayout = Model_t::gridlayout_type; + static constexpr std::size_t dimension = GridLayout::dimension; - return bulkV; + PatchLevel(amr::Hierarchy& hierarchy, Model& model, std::size_t lvl) + : Super(hierarchy, model, lvl) + { } - auto getPopFlux() + auto getB(std::string const& component) { - using Inner = decltype(getPopDensities()); - - std::unordered_map pop_data; - - auto& ions = model_.state.ions; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - for (auto const& pop : ions) - getVecFields(pop.flux(), pop_data, grid, patchID, pop.name()); - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); - - return pop_data; + return this->getVecFieldComponent(this->model.state.electromag.B, + this->vecfield_component(component)); } - auto getEM() + auto getE(std::string const& component) { - using Inner = decltype(getPopDensities()); - - std::unordered_map em_data; - - auto& em = model_.state.electromag; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - for (auto& vecFieldPtr : {&em.B, &em.E}) - { - getVecFields(*vecFieldPtr, em_data, grid, patchID, vecFieldPtr->name()); - } - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, em); - - return em_data; + return this->getVecFieldComponent(this->model.state.electromag.E, + this->vecfield_component(component)); } - - auto getB(std::string componentName) + auto& getIonPop(auto const& popName) { - std::vector, dimension>> patchDatas; - - auto& B = model_.state.electromag.B; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - auto compo = PHARE::core::Components::componentMap().at(componentName); - setPatchDataFromField(patchDatas.emplace_back(), B.getComponent(compo), grid, patchID); - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, B); - - return patchDatas; + for (auto& pop : this->model.state.ions) + if (popName == pop.name()) + return pop; + throw std::runtime_error("no population found named: " + popName); } + auto getNi() { return this->getField(this->model.state.ions.massDensity()); } - auto getE(std::string componentName) + auto getN(std::string const& popName) { - std::vector, dimension>> patchDatas; - - auto& E = model_.state.electromag.E; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - auto compo = PHARE::core::Components::componentMap().at(componentName); - setPatchDataFromField(patchDatas.emplace_back(), E.getComponent(compo), grid, patchID); - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, E); - - return patchDatas; + return this->getField(getIonPop(popName).chargeDensity()); } - - - auto getVi(std::string componentName) + auto getVi(std::string const& component) { - std::vector, dimension>> patchDatas; - - auto& V = model_.state.ions.velocity(); - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - auto compo = PHARE::core::Components::componentMap().at(componentName); - setPatchDataFromField(patchDatas.emplace_back(), V.getComponent(compo), grid, patchID); - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, V); - - return patchDatas; + return this->getVecFieldComponent(this->model.state.ions.velocity(), + this->vecfield_component(component)); } + auto getFlux(std::string const& component, std::string const& popName) + { + return this->getVecFieldComponent(getIonPop(popName).flux(), + this->vecfield_component(component)); + } - - auto getPopFluxCompo(std::string component, std::string popName) + auto getParticles(std::string const& popName) { - std::vector, dimension>> patchDatas; + std::vector*, dimension>> patchDatas; - auto& ions = model_.state.ions; + auto& pop = getIonPop(popName); - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - auto compo = PHARE::core::Components::componentMap().at(component); - for (auto const& pop : ions) - if (pop.name() == popName) - setPatchDataFromField(patchDatas.emplace_back(), pop.flux().getComponent(compo), - grid, patchID); + auto visit = [&](auto& grid, auto patchID, auto /*ilvl*/) { + setPatchDataFromGrid(patchDatas.emplace_back(&pop.domainParticles()), grid, patchID); }; - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); + amr::visitLevel( // + *this->hierarchy.getPatchLevel(this->lvl), this->rm, visit, pop); return patchDatas; } +}; +template +class __attribute__((visibility("hidden"))) +PatchLevel>> : AnyPatchLevel +{ + using Super = AnyPatchLevel; +public: + using Model_t = Model; + using GridLayout = Model_t::gridlayout_type; - auto getEx() { return getE("x"); } - auto getEy() { return getE("y"); } - auto getEz() { return getE("z"); } - - auto getBx() { return getB("x"); } - auto getBy() { return getB("y"); } - auto getBz() { return getB("z"); } - - auto getVix() { return getVi("x"); } - auto getViy() { return getVi("y"); } - auto getViz() { return getVi("z"); } - - auto getFx(std::string pop) { return getPopFluxCompo("x", pop); } - auto getFy(std::string pop) { return getPopFluxCompo("y", pop); } - auto getFz(std::string pop) { return getPopFluxCompo("z", pop); } - - - - - auto getParticles(std::string userPopName) + PatchLevel(amr::Hierarchy& hierarchy, Model& model, std::size_t const lvl) + : Super(hierarchy, model, lvl) { - using Nested = std::vector, dimension>>; - using Inner = std::unordered_map; - - std::unordered_map pop_particles; - - auto getParticleData = [&](Inner& inner, GridLayout& grid, std::string patchID, - std::string key, auto& particles) { - if (particles.size() == 0) - return; - - if (!inner.count(key)) - inner.emplace(key, Nested()); - - auto& patch_data = inner[key].emplace_back(particles.size()); - setPatchDataFromGrid(patch_data, grid, patchID); - core::ParticlePacker{particles}.pack(patch_data.data); - }; - - auto& ions = model_.state.ions; - - auto visit = [&](GridLayout& grid, std::string patchID, std::size_t /*iLevel*/) { - for (auto& pop : ions) - { - if ((userPopName != "" and userPopName == pop.name()) or userPopName == "all") - { - if (!pop_particles.count(pop.name())) - pop_particles.emplace(pop.name(), Inner()); - - auto& inner = pop_particles.at(pop.name()); - - getParticleData(inner, grid, patchID, "domain", pop.domainParticles()); - getParticleData(inner, grid, patchID, "levelGhost", pop.levelGhostParticles()); - } - } - }; - - PHARE::amr::visitLevel(*hierarchy_.getPatchLevel(lvl_), - *model_.resourcesManager, visit, ions); - - return pop_particles; } - -private: - std::size_t lvl_; - amr::Hierarchy& hierarchy_; - HybridModel& model_; }; + } // namespace PHARE::pydata #endif /*PHARE_PYTHON_PATCH_LEVEL_H*/ diff --git a/src/python3/pybind_def.hpp b/src/python3/pybind_def.hpp index a5a9b21ca..614785aa7 100644 --- a/src/python3/pybind_def.hpp +++ b/src/python3/pybind_def.hpp @@ -7,6 +7,7 @@ #include #include "core/utilities/span.hpp" +#include "core/utilities/types.hpp" #include "pybind11/stl.h" #include "pybind11/numpy.h" @@ -35,11 +36,55 @@ std::size_t ndSize(PyArrayInfo const& ar_info) } +template +auto shape(pybind11::buffer_info const& info) +{ + if (info.ndim != dim) + throw std::runtime_error("bad dim"); + return core::for_N_make_array([&](auto i) -> std::size_t { return info.shape[i]; }); +} + +template +auto shape(py_array_t const& ar) +{ + return shape(ar.request()); +} + + +template +auto strides(pybind11::buffer_info const& info) +{ + if (info.ndim != dim) + throw std::runtime_error("bad dim"); + return core::for_N_make_array([&](auto i) -> std::size_t { return info.strides[i]; }); +} + +template +auto strides(py_array_t const& ar) +{ + return strides(ar.request()); +} + +template +std::array strides_from(std::array const& shape) +{ + std::size_t constexpr data_size = sizeof(T); + + if constexpr (dim == 1) + return {data_size}; + + if constexpr (dim == 2) + return {data_size * shape[1], data_size}; + + if constexpr (dim == 3) + return {data_size * shape[1] * shape[2], data_size * shape[1], data_size}; +} + template class __attribute__((visibility("hidden"))) PyArrayWrapper : public core::Span { public: - PyArrayWrapper(PHARE::pydata::py_array_t const& array) + PyArrayWrapper(py_array_t const& array) : core::Span{static_cast(array.request().ptr), pydata::ndSize(array.request())} , _array{array} { @@ -48,7 +93,7 @@ class __attribute__((visibility("hidden"))) PyArrayWrapper : public core::Span _array; + py_array_t _array; }; template @@ -66,6 +111,15 @@ core::Span makeSpan(py_array_t const& py_array) } +template +auto field_as_memory_view(Field& field) +{ + using value_type = Field::value_type; + return pybind11::memoryview::from_buffer(field.data(), field.shape(), + strides_from(field.shape())); +} + + } // namespace PHARE::pydata #endif /*PHARE_PYTHON_PYBIND_DEF_H*/ diff --git a/tests/simulator/CMakeLists.txt b/tests/simulator/CMakeLists.txt index b1226b604..1e275478b 100644 --- a/tests/simulator/CMakeLists.txt +++ b/tests/simulator/CMakeLists.txt @@ -9,8 +9,7 @@ if(NOT ${PHARE_PROJECT_DIR} STREQUAL ${CMAKE_BINARY_DIR}) endif() phare_python3_exec(9 validation test_validation.py ${CMAKE_CURRENT_BINARY_DIR}) -phare_python3_exec(9 data-wrangler data_wrangler.py ${CMAKE_CURRENT_BINARY_DIR}) -phare_python3_exec(9 sim-refineParticlNbr refined_particle_nbr.py ${CMAKE_CURRENT_BINARY_DIR}) +phare_python3_exec(9 data-wrangler test_data_wrangler.py ${CMAKE_CURRENT_BINARY_DIR}) add_no_mpi_python3_test(periodicity test_init_periodicity.py ${CMAKE_CURRENT_BINARY_DIR}) add_no_mpi_python3_test(exceptions test_exceptions.py ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/tests/simulator/refined_particle_nbr.py b/tests/simulator/refined_particle_nbr.py deleted file mode 100644 index 59b3271d3..000000000 --- a/tests/simulator/refined_particle_nbr.py +++ /dev/null @@ -1,132 +0,0 @@ -#!/usr/bin/env python3 -# -# formatted with black - - -import os -import sys -import yaml -import unittest -import numpy as np - -from pyphare import cpp -from pyphare.simulator.simulator import Simulator - -from tests.simulator import NoOverwriteDict -from tests.simulator import populate_simulation -from tests.simulator.config import project_root - - -class SimulatorRefinedParticleNbr(unittest.TestCase): - def __init__(self, *args, **kwargs): - super(SimulatorRefinedParticleNbr, self).__init__(*args, **kwargs) - self.simulator = None - with open(os.path.join(project_root, "res/amr/splitting.yml"), "r") as stream: - try: - self.yaml_root = yaml.safe_load(stream) - except yaml.YAMLError as exc: - print(exc) - sys.exit(1) - - # while splitting particles may leave the domain area - # so we remove the particles from the border cells of each patch - def _less_per_dim(self, dim, refined_particle_nbr, patch): - if dim == 1: - return refined_particle_nbr * 2 - cellNbr = patch.upper - patch.lower + 1 - if dim == 2: - return refined_particle_nbr * ((cellNbr[0] * 2 + (cellNbr[1] * 2))) - raise ValueError("Unhandled dimension for function") - - def _check_deltas_and_weights(self, sim, dim, interp, refined_particle_nbr): - yaml_dim = self.yaml_root["dimension_" + str(dim)] - yaml_interp = yaml_dim["interp_" + str(interp)] - yaml_n_particles = yaml_interp["N_particles_" + str(refined_particle_nbr)] - - yaml_delta = [float(s) for s in str(yaml_n_particles["delta"]).split(" ")] - yaml_weight = [float(s) for s in str(yaml_n_particles["weight"]).split(" ")] - - splitter_t = cpp.splitter_type(sim) - np.testing.assert_allclose(yaml_delta, splitter_t.delta) - np.testing.assert_allclose(yaml_weight, splitter_t.weight) - - def _do_dim(self, dim, min_diff, max_diff): - from pyphare.pharein.simulation import valid_refined_particle_nbr - - for interp in range(1, 4): - prev_split_particle_max = 0 - for refined_particle_nbr in valid_refined_particle_nbr[dim][interp]: - - simInput = NoOverwriteDict( - {"refined_particle_nbr": refined_particle_nbr} - ) - sim = populate_simulation(dim, interp, **simInput) - self._check_deltas_and_weights(sim, dim, interp, refined_particle_nbr) - - self.simulator = Simulator(sim) - self.simulator.initialize() - dw = self.simulator.data_wrangler() - max_per_pop = 0 - leaving_particles = 0 - for pop, particles in dw.getPatchLevel(1).getParticles().items(): - per_pop = 0 - for key, patches in particles.items(): - for patch in patches: - leaving_particles += self._less_per_dim( - dim, refined_particle_nbr, patch - ) - per_pop += patch.data.size() - max_per_pop = max(max_per_pop, per_pop) - prev_min_diff = prev_split_particle_max * min_diff - - # while splitting particles may leave the domain area - # so we remove the particles from the border cells of each patch - self.assertTrue(max_per_pop > prev_min_diff - (leaving_particles)) - if prev_split_particle_max > 0: - prev_max_diff = prev_min_diff * dim * max_diff - self.assertTrue(max_per_pop < prev_max_diff) - prev_split_particle_max = max_per_pop - self.simulator = None - - """ 1d - refine 10 cells in 1d, ppc 100 - 10 * 2 * ppc = 200 - 10 * 3 * ppc = 300 300 / 200 = 1.5 - 10 * 4 * ppc = 400 500 / 400 = 1.33 - 10 * 5 * ppc = 500 500 / 400 = 1.25 - taking the minimul diff across permutations - current to previous should be at least this times bigger - """ - PREVIOUS_ITERATION_MIN_DIFF_1d = 1.25 - PREVIOUS_ITERATION_MAX_DIFF_1d = 1.50 - - def test_1d(self): - This = type(self) - self._do_dim( - 1, This.PREVIOUS_ITERATION_MIN_DIFF_1d, This.PREVIOUS_ITERATION_MAX_DIFF_1d - ) - - """ 2d - refine 10x10 cells in 2d, ppc 100 - 10 * 10 * 4 * ppc = 400 - 10 * 10 * 8 * ppc = 800 800 / 400 = 1.5 - 10 * 10 * 9 * ppc = 900 900 / 800 = 1.125 - """ - PREVIOUS_ITERATION_MIN_DIFF_2d = 1.125 - PREVIOUS_ITERATION_MAX_DIFF_2d = 1.50 - - def test_2d(self): - This = type(self) - self._do_dim( - 2, This.PREVIOUS_ITERATION_MIN_DIFF_2d, This.PREVIOUS_ITERATION_MAX_DIFF_2d - ) - - def tearDown(self): - # needed in case exception is raised in test and Simulator - # not reset properly - if self.simulator is not None: - self.simulator.reset() - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/simulator/test_data_wrangler.py b/tests/simulator/test_data_wrangler.py new file mode 100644 index 000000000..5db83f226 --- /dev/null +++ b/tests/simulator/test_data_wrangler.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 + +import unittest +import numpy as np + +from pyphare import cpp +from pyphare.simulator.simulator import Simulator +from tests.simulator import populate_simulation + + +from pyphare.pharesee.hierarchy.fromsim import hierarchy_from_sim + + +class DataWranglerTest(unittest.TestCase): + def __init__(self, *args, **kwargs): + super(DataWranglerTest, self).__init__(*args, **kwargs) + self.dw = None + self.simulator = None + + def test(self): + ndim, interp = 2, 1 + self.simulator = Simulator(populate_simulation(ndim, interp)) + self.simulator.initialize() + + hier = None + hier = hierarchy_from_sim(self.simulator, qty="EM_B_x", hier=hier, sync=True) + hier = hierarchy_from_sim(self.simulator, qty="EM_B_y", hier=hier, sync=True) + hier = hierarchy_from_sim(self.simulator, qty="density", hier=hier, sync=True) + + if cpp.mpi_rank() == 0: + hier.plot( + filename="data_wrangler.Bx.png", + qty="EM_B_x", + plot_patches=True, + levels=(0,), + ) + hier.plot( + filename="data_wrangler.By.png", + qty="EM_B_y", + plot_patches=True, + levels=(0,), + ) + # hier.plot( # fails? + # filename="data_wrangler.Ni.png", + # qty="density", + # plot_patches=True, + # levels=(0,), + # ) + + def tearDown(self): + del self.dw + if self.simulator is not None: + self.simulator.reset() + + +if __name__ == "__main__": + unittest.main()