divine debugging

src/amr/CMakeLists.txt

@@ -78,6 +78,7 @@ set( SOURCES_CPP
      resources_manager/amr_utils.cpp
      data/field/coarsening/field_coarsen.cpp
      messengers/messenger_factory.cpp
+     debugod.cpp
    )
 
 add_library(${PROJECT_NAME} ${SOURCES_INC} ${SOURCES_CPP})

src/amr/debugod.cpp

@@ -0,0 +1,12 @@
+#include "debugod.hpp"
+
+
+namespace PHARE::amr
+{
+template<typename PHARE_TYPES>
+DEBUGOD<PHARE_TYPES>& DEBUGOD<PHARE_TYPES>::INSTANCE()
+{
+    static DEBUGOD instance;
+    return instance;
+}
+}; // namespace PHARE::amr

src/amr/debugod.hpp

@@ -0,0 +1,216 @@
+#ifndef PHARE_DEBUGOD_HPP
+#define PHARE_DEBUGOD_HPP
+
+#include "core/def.hpp"
+#include "core/utilities/box/box.hpp"
+#include "core/utilities/point/point.hpp"
+#include "core/utilities/mpi_utils.hpp"
+#include "amr/wrappers/hierarchy.hpp"
+#include "amr/data/field/field_data.hpp"
+#include "amr/resources_manager/amr_utils.hpp"
+
+#include <SAMRAI/hier/PatchHierarchy.h>
+
+#include <memory>
+#include <string>
+#include <vector>
+#include <unordered_map>
+#include <cstdint>
+
+namespace PHARE::amr
+{
+
+
+template<typename PHARE_TYPES>
+class DEBUGOD
+{
+public:
+    static constexpr auto dimension    = PHARE_TYPES::dimension;
+    static constexpr auto interp_order = PHARE_TYPES::interp_order;
+    using Grid_t                       = PHARE_TYPES::Grid_t;
+    using GridLayout_t                 = PHARE_TYPES::GridLayout_t;
+    using FieldData_t                  = PHARE::amr::FieldData<GridLayout_t, Grid_t>;
+
+    using Point_t = PHARE::core::Point<double, dimension>;
+
+    struct GodValue
+    {
+        Point_t coords;
+        std::array<int, dimension> loc_index;
+        std::array<int, dimension> amr_index;
+        double value;
+        int rank;
+        int patchID;
+
+        // Add other necessary fields and methods as needed
+    };
+    using GodExtract = std::unordered_map<std::uint32_t, std::vector<GodValue>>;
+
+    NO_DISCARD static DEBUGOD& INSTANCE();
+
+
+    void setHierarchy(std::shared_ptr<SAMRAI::hier::PatchHierarchy> const& hier)
+    {
+        hierarchy_ = hier;
+    }
+
+    bool isActive() const { return hierarchy_ != nullptr; }
+
+    NO_DISCARD auto time_is(std::string name, double time) const
+    {
+        // take first patch of first level
+        // it should be enought to get time
+        // this limits to looking at times at coarser time steps for now
+        auto patch = *(hierarchy_->getPatchLevel(0)->begin());
+        auto pdata = getPatchData(*patch, name);
+        return time == pdata->getTime();
+    }
+
+    NO_DISCARD auto inspect(std::string name, Point_t const& lower, Point_t const& upper) const
+    {
+        GodExtract god_values;
+        for (auto ilvl = 0u; ilvl < hierarchy_->getNumberOfLevels(); ++ilvl)
+        {
+            auto level       = hierarchy_->getPatchLevel(ilvl);
+            god_values[ilvl] = std::vector<GodValue>{};
+
+            for (auto& patch : *level)
+            {
+                if (!is_local(*patch))
+                    continue;
+
+                auto extract_box = PHARE::core::Box<double, dimension>{lower, upper};
+                auto patch_ghost_box
+                    = phare_box_from<dimension, double>(getPatchData(*patch, name)->getGhostBox());
+
+                auto& field    = getField(*patch, name);
+                auto layout    = layoutFromPatch<GridLayout_t>(*patch);
+                auto centering = GridLayout_t::centering(field.physicalQuantity());
+
+                for (auto i = 0u; i < dimension; ++i)
+                {
+                    if (centering[i] == PHARE::core::QtyCentering::primal)
+                    {
+                        extract_box.upper[i] += 1;
+                    }
+                }
+
+                auto intersected_box = patch_ghost_box * extract_box;
+
+                if (!intersected_box)
+                    continue;
+
+
+
+                // loop on nodes
+                // given the mesh_size_ on root level
+                // it is easy to get the level mesh size
+                // and given the lower/upper bounds of the coordinates
+                // it's easy to iterate over all nodes
+                // these if constexpr may be removable
+                // with the FieldBox object maybe....
+
+                GodValue gval;
+                auto box = *intersected_box;
+
+                if constexpr (dimension == 1)
+                {
+                    //
+                }
+
+                else if constexpr (dimension == 2)
+                {
+                    auto& dl     = layout.meshSize();
+                    auto ixStart = static_cast<int>((box.lower[0] - layout.origin()[0]) / dl[0]);
+                    auto ixEnd   = static_cast<int>((box.upper[0] - layout.origin()[0]) / dl[0]);
+                    auto iyStart = static_cast<int>((box.lower[1] - layout.origin()[1]) / dl[1]);
+                    auto iyEnd   = static_cast<int>((box.upper[1] - layout.origin()[1]) / dl[1]);
+
+                    for (auto ix = ixStart; ix <= ixEnd; ++ix)
+                    {
+                        for (auto iy = iyStart; iy <= iyEnd; ++iy)
+                        {
+                            gval.coords
+                                = layout.fieldNodeCoordinates(field, layout.origin(), ix, iy);
+                            gval.value = field(ix, iy);
+                        }
+                    }
+                }
+                else if constexpr (dimension == 3)
+                {
+                    // for (auto& node : intersected_box)
+                    // {
+                    // }
+                }
+                god_values[ilvl].push_back(gval);
+            }
+        }
+
+        return god_values;
+    }
+
+
+    NO_DISCARD auto inspect(std::string name, Point_t const& coord)
+    {
+        return inspect(name, coord, coord);
+    }
+
+
+
+    void print(GodExtract const& god_values)
+    {
+        for (auto& [ilvl, values] : god_values)
+        {
+            std::cout << "Level " << ilvl << ":\n";
+            for (auto& v : values)
+            {
+                auto& coords  = v.coords;
+                auto& loc_idx = v.loc_index;
+                auto& amr_idx = v.loc_index;
+                auto& rank    = v.rank;
+                auto& patchID = v.patchID;
+
+                std::cout << "\n";
+            }
+        }
+    }
+
+
+    // void stop() { god_.release(); }
+
+    // NO_DISCARD auto& god()
+    // {
+    //     if (!god_)
+    //         init();
+    //     return *god_;
+    // }
+
+private:
+    bool is_local(SAMRAI::hier::Patch const& patch) const
+    {
+        return patch.getBox().getBoxId().getOwnerRank() == PHARE::core::mpi::rank();
+    }
+
+    auto getPatchData(SAMRAI::hier::Patch const& patch, std::string const& name) const
+    {
+        auto db      = SAMRAI::hier::VariableDatabase::getDatabase();
+        auto var_id  = db->getVariable(name);
+        auto context = db->getContext("default");
+        return patch.getPatchData(var_id, context);
+    }
+
+    auto& getField(SAMRAI::hier::Patch const& patch, std::string const& name) const
+    {
+        auto pdata            = getPatchData(patch, name);
+        auto const& fielddata = std::dynamic_pointer_cast<FieldData_t>(pdata);
+        return fielddata->field;
+    }
+
+
+
+    DEBUGOD() {}
+    std::shared_ptr<SAMRAI::hier::PatchHierarchy> hierarchy_;
+};
+}; // namespace PHARE::amr
+
+#endif

src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp

@@ -24,6 +24,9 @@
 #include "amr/resources_manager/amr_utils.hpp"
 #include "amr/data/field/refine/magnetic_refine_patch_strategy.hpp"
 
+#include "amr/debugod.hpp"
+#include "phare_core.hpp"
+
 #include "core/numerics/interpolator/interpolator.hpp"
 #include "core/hybrid/hybrid_quantities.hpp"
 #include "core/data/particles/particle_array.hpp"
@@ -647,6 +650,17 @@ namespace amr
             magSharedNodeRefineSchedules[levelNumber]->fillData(time);
             elecSharedNodesRefiners_.fill(hybridModel.state.electromag.E, levelNumber, time);
 
+
+            auto& god = DEBUGOD<PHARE::core::PHARE_Types<2, 1>>::INSTANCE();
+            if (god.isActive())
+            {
+                if (god.time_is("EM_B_x", 0.225))
+                {
+                    auto bx_dbg = god.inspect("EM_B_x", {12.2, 13.5});
+                    god.print(bx_dbg);
+                    auto bx_dbg_rge = god.inspect("EM_B_x", {12.2, 8.0}, {12.6, 9.});
+                }
+            }
             // we fill magnetic field ghosts only on patch ghost nodes and not on level
             // ghosts the reason is that 1/ filling ghosts is necessary to prevent mismatch
             // between ghost and overlaped neighboring patch domain nodes resulting from

src/amr/utilities/box/amr_box.hpp

@@ -24,10 +24,12 @@ NO_DISCARD auto samrai_box_from(PHARE::core::Box<Type, dim> const& box, int samr
 template<std::size_t dim, typename Type = int>
 NO_DISCARD auto phare_box_from(SAMRAI::hier::Box const& box)
 {
-    std::array<Type, dim> lower = *reinterpret_cast<std::array<int, dim> const*>(&box.lower()[0]);
-    std::array<Type, dim> upper = *reinterpret_cast<std::array<int, dim> const*>(&box.upper()[0]);
+    std::array<int, dim> lower = *reinterpret_cast<std::array<int, dim> const*>(&box.lower()[0]);
+    std::array<int, dim> upper = *reinterpret_cast<std::array<int, dim> const*>(&box.upper()[0]);
 
-    return PHARE::core::Box<Type, dim>{core::Point{lower}, core::Point{upper}};
+    auto gen = [](auto v) { return static_cast<Type>(v); };
+    return PHARE::core::Box<Type, dim>{core::Point{core::generate(gen, lower)},
+                                       core::Point{core::generate(gen, upper)}};
 }
 
 NO_DISCARD inline bool operator==(SAMRAI::hier::Box const& b1, SAMRAI::hier::Box const& b2)

src/core/utilities/box/box.hpp

@@ -25,7 +25,7 @@ class box_iterator;
 template<typename Type, std::size_t dim>
 struct Box
 {
-    static const size_t dimension = dim;
+    static size_t const dimension = dim;
 
 
     Point<Type, dim> lower;

src/simulator/simulator.hpp

@@ -16,6 +16,7 @@
 #include "amr/load_balancing/load_balancer_details.hpp"
 #include "amr/load_balancing/load_balancer_manager.hpp"
 #include "amr/load_balancing/load_balancer_estimator_hybrid.hpp"
+#include "amr/debugod.hpp"
 
 namespace PHARE
 {
@@ -339,6 +340,8 @@ Simulator<_dimension, _interp_order, _nbRefinedPart>::Simulator(
     , functors_{functors_setup(dict)}
     , multiphysInteg_{std::make_shared<MultiPhysicsIntegrator>(dict["simulation"], functors_)}
 {
+    PHARE::amr::DEBUGOD<core::PHARE_Types<_dimension, _interp_order>>::INSTANCE().setHierarchy(
+        hierarchy_);
     if (find_model("HybridModel"))
         hybrid_init(dict);
     else