more

pyphare/pyphare/core/box.py

@@ -68,6 +68,34 @@ def __sub__(self, other):
     def copy(self):
         return Box(self.lower.copy(), self.upper.copy())
 
+    def __iter__(self):
+        return BoxIterator(self)
+
+
+class BoxIterator:
+    def __init__(self, box):
+        self.l = box.lower
+        self.u = box.upper
+        self.c = l.copy()
+        self.c[0] -= 1  # hax for first iteration
+
+    def __iter__(self):
+        return self
+
+    def __next__(self):
+        if (self.c >= self.u).all():
+            raise StopIteration()
+        self.c[0] += 1
+        if len(self.c) > 1:
+            if self.c[0] > self.u[0]:
+                self.c[1] += 1
+                self.c[0] = self.l[0]
+        if len(self.c) > 2:
+            if self.c[1] > self.u[1]:
+                self.c[2] += 1
+                self.c[1] = self.l[1]
+        return self.c
+
 
 class nDBox(Box):
     def __init__(self, dim, l, u):

pyphare/pyphare/pharesee/particles.py

@@ -75,6 +75,14 @@ def shift_icell(self, offset):
     def size(self):
         return len(self.weights)
 
+    def __ne__(self, that):
+        assert isinstance(that, Particles)
+        try:
+            return not self == that
+        except AssertionError:
+            return True
+
+
     def __eq__(self, that):
         if isinstance(that, Particles):
             # fails on OSX for some reason

src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp

@@ -70,6 +70,8 @@ void ConcreteLoadBalancerHybridStrategyHomogeneous<PHARE_T>::compute(
         // of ghost cell for this patch data is null.
 
         // The lb_view is a CellData, meaning that it is dual as the index of an amr box
+
+        // todo replace below with std::fill/std::fill_n data, size = 1
         core::Box<std::uint32_t, dimension> local_box{
             core::Point{core::ConstArray<std::uint32_t, dimension>()},
             core::Point{

src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp

@@ -109,7 +109,7 @@ namespace amr
         using DefaultCoarsenOp  = BaseCoarsenOp<DefaultFieldCoarsener<dimension>>;
 
     public:
-        static const std::string stratName;
+        static const inline std::string stratName    = "HybridModel-HybridModel";
         static constexpr std::size_t rootLevelNumber = 0;
 
 
@@ -1086,9 +1086,6 @@ namespace amr
         CoarsenOperator_ptr magneticCoarseningOp_{std::make_shared<MagneticCoarsenOp>()};
     };
 
-    template<typename HybridModel, typename RefinementParams>
-    const std::string HybridHybridMessengerStrategy<HybridModel, RefinementParams>::stratName
-        = "HybridModel-HybridModel";
 
 } // namespace amr
 

src/amr/multiphysics_integrator.hpp

@@ -30,6 +30,7 @@
 #include "amr/solvers/solver_mhd.hpp"
 #include "amr/solvers/solver_ppc.hpp"
 
+#include "core/logger.hpp"
 #include "core/utilities/algorithm.hpp"
 
 #include "load_balancing/load_balancer_manager.hpp"
@@ -300,7 +301,7 @@ namespace solver
          */
         void initializeLevelData(std::shared_ptr<SAMRAI::hier::PatchHierarchy> const& hierarchy,
                                  int const levelNumber, double const initDataTime,
-                                 bool const /*canBeRefined*/, bool const /*initialTime*/,
+                                 bool const canBeRefined, bool const initialTime,
                                  std::shared_ptr<SAMRAI::hier::PatchLevel> const& oldLevel
                                  = std::shared_ptr<SAMRAI::hier::PatchLevel>(),
                                  bool const allocateData = true) override
@@ -311,14 +312,18 @@ namespace solver
             auto& levelInitializer = getLevelInitializer(model.name());
 
             bool const isRegridding = oldLevel != nullptr;
-            auto level              = hierarchy->getPatchLevel(levelNumber);
+            // auto level              = hierarchy->getPatchLevel(levelNumber);
+
+            PHARE_LOG_LINE_STR("init level " << levelNumber //
+                                             << " with regriding = " << isRegridding
+                                             << " with allocateData = " << allocateData
+                                             << " with canBeRefined = " << canBeRefined
+                                             << " with initialTime = " << initialTime);
 
-            std::cout << "init level " << levelNumber << " with regriding = " << isRegridding
-                      << "\n";
             PHARE_LOG_START("initializeLevelData::allocate block");
             if (allocateData)
             {
-                for (auto patch : *level)
+                for (auto patch : *hierarchy->getPatchLevel(levelNumber))
                 {
                     model.allocate(*patch, initDataTime);
                     solver.allocate(model, *patch, initDataTime);

src/amr/wrappers/integrator.hpp

@@ -20,7 +20,6 @@
 #include <SAMRAI/tbox/MemoryDatabase.h>
 
 
-
 #include "initializer/data_provider.hpp"
 
 
@@ -82,9 +81,6 @@ Integrator<_dimension>::Integrator(
     double endTime)
     : rebalance_coarsest{check_rebalance_coarsest(dict)}
 {
-    // auto loadBalancer = std::make_shared<SAMRAI::mesh::TreeLoadBalancer>(
-    //     SAMRAI::tbox::Dimension{dimension}, "LoadBalancer");
-
     auto loadBalancer_db = std::make_shared<SAMRAI::tbox::MemoryDatabase>("LoadBalancerDB");
     loadBalancer_db->putDouble("flexible_load_tolerance", .75);
     auto loadBalancer = std::make_shared<SAMRAI::mesh::CascadePartitioner>(

tests/simulator/test_load_balancing.py

@@ -4,9 +4,12 @@
 #   historgrams?
 #
 
+import pyphare.core.box as boxm
 import pyphare.pharein as ph
 from pyphare.simulator.simulator import Simulator, startMPI
 
+from pyphare.pharesee.particles import aggregate as merge_particles
+
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib as mpl
@@ -18,12 +21,14 @@
 cpp = cpp_lib()
 startMPI()
 
-time_step_nbr = 1000
-time_step = 0.001
+time_step_nbr = 1
+time_step = 0.005
 smallest_patch_size = 10
 largest_patch_size = 20
 rebalance_coarsest = False
+cells = (100, 100)
 diag_outputs = "phare_outputs/harris/2d/load_balancing"
+timestamps = [0, 0.005]
 
 
 def config():
@@ -32,7 +37,7 @@ def config():
         largest_patch_size=largest_patch_size,
         time_step_nbr=time_step_nbr,
         time_step=time_step,
-        cells=(100, 100),
+        cells=cells,
         dl=(0.2, 0.2),
         refinement_boxes={},
         hyper_resistivity=0.001,
@@ -137,7 +142,7 @@ def vthz(x, y):
     ph.ElectronModel(closure="isothermal", Te=0.0)
     ph.ParticleDiagnostics(
         quantity="domain",
-        write_timestamps=[0, sim.final_time],
+        write_timestamps=timestamps,
         population_name="protons",
     )
     return sim
@@ -152,36 +157,71 @@ def get_time(path, time, datahier=None):
     )
 
 
-def post_advance(new_time):
-    if cpp.mpi_rank() == 0:
-        print(f"running tests at time {new_time}")
-        from tests.simulator.test_advance import AdvanceTestBase
+def check_time(sim, time=0):
+    print("check_time", time)
+    hier = get_time(diag_outputs, time)
+    assert len(hier.levels()) == 1  # L0 ONLY!
 
-        test = AdvanceTestBase()
-        test.base_test_overlaped_fields_are_equal(
-            get_time(diag_outputs, new_time), new_time
-        )
-        print(f"tests passed")
+    ppc = np.ones(np.asarray(sim.cells))
 
+    for ilvl, lvl in hier.levels().items():
+        for patch in lvl:
+            for pd_key, pd in patch.patch_datas.items():
+                for bit in bt(patch.box):
+                    box = boxm.Box(bit, bit)
+                    slice = pd.dataset.select(box)
+                    ppc[tuple(bit)] = slice.size()
+    return hier, ppc
 
-def check_time(sim, time=0):
-    print("check_time", time)
+
+def get_merged_particles(sim, time=0):
+    print("get_merged_particles", time)
     hier = get_time(diag_outputs, time)
     assert len(hier.levels()) == 1  # L0 ONLY!
+
+    particles = []
     for ilvl, lvl in hier.levels().items():
         for patch in lvl:
             for pd_key, pd in patch.patch_datas.items():
-                print("len(pd.dataset)", pd.dataset.size())
+                particles += [pd.dataset]
+
+    assert len(particles) > 1
+    return merge_particles(particles)
+
+
+def make_fig(hier, fig_name, ilvl, collections):
+    hier.plot_2d_patches(0, collections=collections).savefig(fig_name + ".png")
 
 
-def test_balance(sim):
-    check_time(sim)
-    check_time(sim, time_step_nbr * time_step)
+def test_balance(sim, time):
+    hier, ppc = check_time(sim, time)
+    max_ppc = np.max(ppc)
+    ppc_collections = [
+        {
+            "boxes": [boxm.Box(index, index)],
+            "facecolor": "purple",
+            "alpha": value / max_ppc,
+        }
+        for index, value in np.ndenumerate(ppc)
+    ]
+    make_fig(hier, f"lb_t{time}", 0, ppc_collections)
+
+
+def test_particles_have_evolved(sim):
+    init = get_merged_particles(sim)
+    end = get_merged_particles(sim, time_step)
+
+    assert init != end
 
 
 def main():
-    sim = Simulator(config()).run()
-    test_balance(sim)
+    sim = config()
+    Simulator(sim).run()
+
+    test_particles_have_evolved(sim)
+
+    # test_balance(sim, timestamps[cpp.mpi_rank() + 1])
+    # test_balance(sim, 0)
 
 
 if __name__ == "__main__":

tests/simulator/test_restarts.py

@@ -124,6 +124,7 @@ def vthxyz(x):
     dl=0.3,
     diag_options=dict(format="phareh5", options=dict(dir=out, mode="overwrite")),
     restart_options=dict(dir=out, mode="overwrite"),
+    advanced={"integrator/rebalance_coarsest": False},
 )