advanced simulation paramaterization

pyphare/pyphare/pharein/__init__.py

@@ -215,6 +215,13 @@ def as_paths(rb):
     add_double("simulation/algo/ohm/resistivity", simulation.resistivity)
     add_double("simulation/algo/ohm/hyper_resistivity", simulation.hyper_resistivity)
 
+    for k, v in simulation.advanced.items():
+        path = f"simulation/advanced/{k}"
+        if isinstance(v, int):
+            add_int(path, v)
+        else:
+            add_string(path, v)
+
     init_model = simulation.model
     modelDict = init_model.model_dict
 

pyphare/pyphare/pharein/simulation.py

@@ -616,6 +616,7 @@ def wrapper(simulation_object, **kwargs):
             "tag_buffer",
             "description",
             "loadbalancing",
+            "advanced",
         ]
 
         accepted_keywords += check_optional_keywords(**kwargs)
@@ -684,6 +685,8 @@ def wrapper(simulation_object, **kwargs):
 
         kwargs["hyper_resistivity"] = check_hyper_resistivity(**kwargs)
 
+        kwargs["advanced"] = kwargs.get("advanced", {})
+
         return func(simulation_object, **kwargs)
 
     return wrapper

src/amr/wrappers/integrator.hpp

@@ -13,7 +13,6 @@
 #include <SAMRAI/mesh/BergerRigoutsos.h>
 #include <SAMRAI/mesh/GriddingAlgorithm.h>
 #include <SAMRAI/mesh/StandardTagAndInitialize.h>
-// #include <SAMRAI/mesh/TreeLoadBalancer.h>
 #include <SAMRAI/mesh/CascadePartitioner.h>
 #include <SAMRAI/tbox/Database.h>
 #include <SAMRAI/tbox/DatabaseBox.h>
@@ -30,7 +29,12 @@ namespace PHARE::amr
 template<std::size_t _dimension>
 class Integrator
 {
-    bool static constexpr rebalance_coarsest = false;
+    bool check_rebalance_coarsest(initializer::PHAREDict const& dict)
+    {
+        return initializer::dict_get(dict, "simulation/advanced/integrator/rebalance_coarsest",
+                                     int{1})
+               > 0;
+    }
 
 public:
     static constexpr std::size_t dimension = _dimension;
@@ -50,6 +54,7 @@ class Integrator
                double startTime, double endTime);
 
 private:
+    bool rebalance_coarsest;
     std::shared_ptr<SAMRAI::algs::TimeRefinementIntegrator> timeRefIntegrator_;
 };
 
@@ -75,10 +80,13 @@ Integrator<_dimension>::Integrator(
     std::shared_ptr<SAMRAI::algs::TimeRefinementLevelStrategy> timeRefLevelStrategy,
     std::shared_ptr<SAMRAI::mesh::StandardTagAndInitStrategy> tagAndInitStrategy, double startTime,
     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>(
         SAMRAI::tbox::Dimension{dimension}, "LoadBalancer");
 

src/initializer/data_provider.hpp

@@ -8,7 +8,9 @@
 #include <functional>
 
 #include "cppdict/include/dict.hpp"
+
 #include "core/def.hpp"
+#include "core/logger.hpp"
 #include "core/utilities/span.hpp"
 
 namespace PHARE
@@ -69,6 +71,7 @@ namespace initializer
             return *phareDict;
         }
 
+
     private:
         PHAREDictHandler() = default;
         std::unique_ptr<PHAREDict> phareDict;
@@ -89,4 +92,53 @@ namespace initializer
 } // namespace initializer
 } // namespace PHARE
 
+namespace PHARE::initializer
+{
+
+template<typename Dict, typename Paths>
+std::optional<Dict> _traverse_to_node(Dict const& dict, Paths const& paths, std::size_t idx = 0)
+{
+    if (dict.contains(paths[idx]))
+    {
+        auto const& next = dict[paths[idx]];
+        if (idx == paths.size() - 1)
+            return next;
+        return _traverse_to_node(next, paths, ++idx);
+    }
+    return std::nullopt;
+}
+
+template<typename Dict>
+std::optional<PHAREDict> traverse_to_node(Dict const& dict, std::string const& path)
+{
+    char delimiter  = '/';
+    std::string tmp = "";
+    std::vector<std::string> paths;
+    std::istringstream iss(path);
+    while (std::getline(iss, tmp, delimiter))
+        paths.push_back(tmp);
+    return _traverse_to_node(dict, paths);
+}
+
+template<typename T, typename Dict>
+T const& dict_get(Dict const& dict, std::string const& path)
+{
+    // error if not found
+    auto leaf = traverse_to_node(dict, path);
+    if (leaf)
+        return leaf->template to<T>();
+    throw std::runtime_error("Dict contains no path " + path);
+}
+
+template<typename T, typename Dict>
+T const& dict_get(Dict const& dict, std::string const& path, T const& default_value)
+{
+    auto leaf = traverse_to_node(dict, path);
+    if (leaf)
+        return leaf->template to<T>();
+    return default_value;
+}
+
+} // namespace PHARE::initializer
+
 #endif // DATA_PROVIDER_HPP

tests/simulator/test_load_balancing.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+#
+# TODO: finish
+#   historgrams?
+#
+
+import pyphare.pharein as ph
+from pyphare.simulator.simulator import Simulator, startMPI
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+mpl.use("Agg")
+
+from pyphare.cpp import cpp_lib
+
+cpp = cpp_lib()
+startMPI()
+
+time_step_nbr = 1000
+time_step = 0.001
+smallest_patch_size = 10
+largest_patch_size = 20
+rebalance_coarsest = False
+diag_outputs = "phare_outputs/harris/2d/load_balancing"
+
+
+def config():
+    sim = ph.Simulation(
+        smallest_patch_size=smallest_patch_size,
+        largest_patch_size=largest_patch_size,
+        time_step_nbr=time_step_nbr,
+        time_step=time_step,
+        cells=(100, 100),
+        dl=(0.2, 0.2),
+        refinement_boxes={},
+        hyper_resistivity=0.001,
+        resistivity=0.001,
+        diag_options={
+            "format": "phareh5",
+            "options": {"dir": diag_outputs, "mode": "overwrite"},
+        },
+        advanced={"integrator/rebalance_coarsest": rebalance_coarsest},
+    )
+
+    def density(x, y):
+        L = sim.simulation_domain()[1]
+        return (
+            0.2
+            + 1.0 / np.cosh((y - L * 0.3) / 0.5) ** 2
+            + 1.0 / np.cosh((y - L * 0.7) / 0.5) ** 2
+        )
+
+    def S(y, y0, l):
+        return 0.5 * (1.0 + np.tanh((y - y0) / l))
+
+    def by(x, y):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        w1 = 0.2
+        w2 = 1.0
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+        w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2))
+        w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2))
+        w5 = 2.0 * w1 / w2
+        return (w5 * x0 * w3) + (-w5 * x0 * w4)
+
+    def bx(x, y):
+        Lx = sim.simulation_domain()[0]
+        Ly = sim.simulation_domain()[1]
+        w1 = 0.2
+        w2 = 1.0
+        x0 = x - 0.5 * Lx
+        y1 = y - 0.3 * Ly
+        y2 = y - 0.7 * Ly
+        w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2))
+        w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2))
+        w5 = 2.0 * w1 / w2
+        v1 = -1
+        v2 = 1.0
+        return (
+            v1
+            + (v2 - v1) * (S(y, Ly * 0.3, 0.5) - S(y, Ly * 0.7, 0.5))
+            + (-w5 * y1 * w3)
+            + (+w5 * y2 * w4)
+        )
+
+    def bz(x, y):
+        return 0.0
+
+    def b2(x, y):
+        return bx(x, y) ** 2 + by(x, y) ** 2 + bz(x, y) ** 2
+
+    def T(x, y):
+        K = 1
+        temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5)
+        assert np.all(temp > 0)
+        return temp
+
+    def vx(x, y):
+        return 0.0
+
+    def vy(x, y):
+        return 0.0
+
+    def vz(x, y):
+        return 0.0
+
+    def vthx(x, y):
+        return np.sqrt(T(x, y))
+
+    def vthy(x, y):
+        return np.sqrt(T(x, y))
+
+    def vthz(x, y):
+        return np.sqrt(T(x, y))
+
+    vvv = {
+        "vbulkx": vx,
+        "vbulky": vy,
+        "vbulkz": vz,
+        "vthx": vthx,
+        "vthy": vthy,
+        "vthz": vthz,
+        "nbr_part_per_cell": 100,
+    }
+
+    ph.MaxwellianFluidModel(
+        bx=bx,
+        by=by,
+        bz=bz,
+        protons={"charge": 1, "density": density, **vvv, "init": {"seed": 12334}},
+    )
+    ph.ElectronModel(closure="isothermal", Te=0.0)
+    ph.ParticleDiagnostics(
+        quantity="domain",
+        write_timestamps=[0, sim.final_time],
+        population_name="protons",
+    )
+    return sim
+
+
+def get_time(path, time, datahier=None):
+    time = "{:.10f}".format(time)
+    from pyphare.pharesee.hierarchy import hierarchy_from
+
+    return hierarchy_from(
+        h5_filename=path + "/ions_pop_protons_domain.h5", time=time, hier=datahier
+    )
+
+
+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
+
+        test = AdvanceTestBase()
+        test.base_test_overlaped_fields_are_equal(
+            get_time(diag_outputs, new_time), new_time
+        )
+        print(f"tests passed")
+
+
+def check_time(sim, time=0):
+    print("check_time", time)
+    hier = get_time(diag_outputs, time)
+    assert len(hier.levels()) == 1  # L0 ONLY!
+    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())
+
+
+def test_balance(sim):
+    check_time(sim)
+    check_time(sim, time_step_nbr * time_step)
+
+
+def main():
+    sim = Simulator(config()).run()
+    test_balance(sim)
+
+
+if __name__ == "__main__":
+    main()