update restart tests to allow ndim > 1d

src/amr/multiphysics_integrator.hpp

@@ -645,20 +645,6 @@ namespace solver
 
 
 
-        // bool existloadBalancererOnRange_(int coarsestLevel, int finestLevel)
-        // {
-        //     for (auto iLevel = coarsestLevel; iLevel <= finestLevel; ++iLevel)
-        //     {
-        //         if (levelDescriptors_[iLevel].loadBalancerIndex != LevelDescriptor::NOT_SET)
-        //         {
-        //             return true;
-        //         }
-        //     }
-        //     return false;
-        // }
-
-
-
         bool existModelOnRange_(int coarsestLevel, int finestLevel)
         {
             bool hasModel = true;

src/amr/wrappers/integrator.hpp

@@ -29,21 +29,29 @@ namespace PHARE::amr
 template<std::size_t _dimension>
 class Integrator
 {
-    int static constexpr rebalance_coarsest_default               = 0;
-    std::size_t static constexpr rebalance_coarsest_every_default = 1000;
+    int static constexpr rebalance_coarsest_every_default = 1000;
 
     bool static _rebalance_coarsest(initializer::PHAREDict const& dict)
     {
-        return initializer::dict_get(dict, "simulation/advanced/integrator/rebalance_coarsest",
-                                     rebalance_coarsest_default)
+        return initializer::dict_get(dict, "simulation/advanced/integrator/rebalance_coarsest", 0)
                > 0;
     }
 
-    std::size_t static _rebalance_coarsest_every(initializer::PHAREDict const& dict)
+    bool static _rebalance_coarsest_on_init(initializer::PHAREDict const& dict)
     {
         return initializer::dict_get(dict,
-                                     "simulation/advanced/integrator/rebalance_coarsest_every",
-                                     rebalance_coarsest_every_default);
+                                     "simulation/advanced/integrator/rebalance_coarsest_on_init", 0)
+               > 0;
+    }
+
+    std::size_t static _rebalance_coarsest_every(initializer::PHAREDict const& dict)
+    {
+        auto in
+            = initializer::dict_get(dict, "simulation/advanced/integrator/rebalance_coarsest_every",
+                                    rebalance_coarsest_every_default);
+        if (in < 0)
+            throw std::runtime_error("rebalance_coarsest_every must be positive");
+        return static_cast<std::size_t>(in);
     }
 
     bool static _is_tagging_refinement(initializer::PHAREDict const& dict)
@@ -53,8 +61,6 @@ class Integrator
                == std::string{"auto"};
     }
 
-
-
 public:
     static constexpr std::size_t dimension = _dimension;
 
@@ -65,6 +71,11 @@ class Integrator
                                            or (time_step_idx > 0 and rebalance_coarsest_every > 0
                                                and time_step_idx % rebalance_coarsest_every == 0));
 
+        PHARE_LOG_LINE_STR(is_tagging_refinement
+                           << " " << time_step_idx << " " << rebalance_coarsest << " "
+                           << rebalance_coarsest_on_init << " " << rebalance_coarsest_every << " "
+                           << rebalance_coarsest_now);
+
         auto new_time = timeRefIntegrator_->advanceHierarchy(dt, rebalance_coarsest_now);
         ++time_step_idx;
         return new_time;
@@ -81,8 +92,8 @@ class Integrator
 
 private:
     bool is_tagging_refinement                 = false;
-    bool rebalance_coarsest                    = true;
-    bool rebalance_coarsest_on_init            = true;
+    bool rebalance_coarsest                    = false;
+    bool rebalance_coarsest_on_init            = false;
     std::size_t time_step_idx                  = 0;
     std::size_t const rebalance_coarsest_every = rebalance_coarsest_every_default;
 
@@ -115,6 +126,7 @@ Integrator<_dimension>::Integrator(
     double startTime, double endTime)
     : is_tagging_refinement{_is_tagging_refinement(dict)}
     , rebalance_coarsest{_rebalance_coarsest(dict)}
+    , rebalance_coarsest_on_init{_rebalance_coarsest_on_init(dict)}
     , rebalance_coarsest_every{_rebalance_coarsest_every(dict)}
 {
     loadBalancer->setSAMRAI_MPI(

tests/functional/harris/harris_2d.py

@@ -1,7 +1,10 @@
 #!/usr/bin/env python3
 
+from datetime import datetime
+
 import pyphare.pharein as ph
 from pyphare.simulator.simulator import Simulator, startMPI
+from pyphare.pharesee.run import Run
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -14,27 +17,37 @@
 cpp = cpp_lib()
 startMPI()
 
-diag_outputs = "phare_outputs/test/harris/2d"
-from datetime import datetime
+diag_outputs0 = "phare_outputs/balanceL1/harris/2d"
+diag_outputs1 = "phare_outputs/balanceL1/harris/2d_rebal"
+
+time_step = 0.001
+time_step_nbr = 30000
+final_time = time_step * time_step_nbr
 
+# array([10., 20., 30., 40., 50.])
+timestamps = np.arange(
+    0, final_time + time_step, time_step * time_step_nbr / 5, dtype=float
+)[1:]
 
-def config():
+
+def config(diag_dig, **kwargs):
     sim = ph.Simulation(
+        dl=(0.2, 0.2),
+        cells=(100, 100),
+        time_step=time_step,
+        time_step_nbr=time_step_nbr,
         smallest_patch_size=15,
         largest_patch_size=25,
-        time_step_nbr=1000,
-        time_step=0.001,
-        # boundary_types="periodic",
-        cells=(100, 100),
-        dl=(0.2, 0.2),
-        refinement_boxes={},
+        max_nbr_levels=2,
+        refinement="tagging",
         hyper_resistivity=0.001,
         resistivity=0.001,
         diag_options={
             "format": "phareh5",
-            "options": {"dir": diag_outputs, "mode": "overwrite"},
+            "options": {"dir": diag_dig, "mode": "overwrite"},
         },
-        strict=True,
+        restart_options=dict(dir=diag_dig, mode="overwrite", timestamps=timestamps),
+        **kwargs,
     )
 
     def density(x, y):
@@ -93,31 +106,19 @@ def T(x, y):
         assert np.all(temp > 0)
         return temp
 
-    def vx(x, y):
-        return 0.0
-
-    def vy(x, y):
+    def vxyz(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):
+    def vthxyz(x, y):
         return np.sqrt(T(x, y))
 
     vvv = {
-        "vbulkx": vx,
-        "vbulky": vy,
-        "vbulkz": vz,
-        "vthx": vthx,
-        "vthy": vthy,
-        "vthz": vthz,
+        "vbulkx": vxyz,
+        "vbulky": vxyz,
+        "vbulkz": vxyz,
+        "vthx": vthxyz,
+        "vthy": vthxyz,
+        "vthz": vthxyz,
         "nbr_part_per_cell": 100,
     }
 
@@ -127,49 +128,73 @@ def vthz(x, y):
         bz=bz,
         protons={"charge": 1, "density": density, **vvv, "init": {"seed": 12334}},
     )
-
     ph.ElectronModel(closure="isothermal", Te=0.0)
-
-    dt = 10 * sim.time_step
-    nt = sim.final_time / dt + 1
-    timestamps = dt * np.arange(nt)
-
     for quantity in ["E", "B"]:
         ph.ElectromagDiagnostics(
             quantity=quantity,
             write_timestamps=timestamps,
             compute_timestamps=timestamps,
         )
-
+    ph.global_vars.sim = None
     return sim
 
 
-def get_time(path, time, datahier=None):
-    time = "{:.10f}".format(time)
-    from pyphare.pharesee.hierarchy import hierarchy_from
-
-    datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", time=time, hier=datahier)
-    datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", time=time, hier=datahier)
-    return datahier
-
-
-def post_advance(new_time):
+def plot(diag_dir):
     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")
+        run = Run(diag_dir)
+        for time in timestamps:
+            run.GetDivB(time).plot(
+                filename=f"{diag_dir}/harris_divb_t{time}.png", plot_patches=True
+            )
+            run.GetB(time).plot(
+                filename=f"{diag_dir}/harris_bx_t{time}.png",
+                qty="Bx",
+                plot_patches=True,
+            )
+            run.GetB(time).plot(
+                filename=f"{diag_dir}/harris_by_t{time}.png",
+                qty="By",
+                plot_patches=True,
+            )
+            run.GetB(time).plot(
+                filename=f"{diag_dir}/harris_bz_t{time}.png",
+                qty="Bz",
+                plot_patches=True,
+            )
+
+            run.GetJ(time).plot(
+                filename=f"{diag_dir}/harris_Jz_t{time}.png",
+                qty="Jz",
+                plot_patches=True,
+                vmin=-2,
+                vmax=2,
+            )
+
+    cpp.mpi_barrier()
+
+
+def run(sim, diag_dir):
+    ph.global_vars.sim = sim
+
+    Simulator(sim).run().reset()
+
+    plot(diag_dir)
+
+    ph.global_vars.sim = None
 
 
 def main():
-    s = Simulator(config(), post_advance=post_advance)
-    s.initialize()
-    post_advance(0)
-    s.run()
+    sim0 = config(diag_outputs0)
+    sim1 = config(
+        diag_outputs1,
+        advanced={
+            "integrator/rebalance_coarsest": 1,
+            "integrator/rebalance_coarsest_every": 5000,
+            "integrator/flexible_load_tolerance": 0.5,
+        },
+    )
+    # run(sim0, diag_outputs0)
+    run(sim1, diag_outputs1)
 
 
 if __name__ == "__main__":

tests/simulator/__init__.py

@@ -1,7 +1,10 @@
 import unittest
+
+from copy import deepcopy
+
+
 from datetime import datetime
 import pyphare.pharein as ph, numpy as np
-from pyphare.pharein import ElectronModel
 
 
 def parse_cli_args(pop_from_sys=True):
@@ -30,32 +33,35 @@ def basicSimulatorArgs(dim: int, interp: int, **kwargs):
     from pyphare.pharein.simulation import valid_refined_particle_nbr
     from pyphare.pharein.simulation import check_patch_size
 
+    args = deepcopy(kwargs)
+
     cells = kwargs.get("cells", [20 for i in range(dim)])
     if not isinstance(cells, (list, tuple)):
         cells = [cells] * dim
 
     _, smallest_patch_size = check_patch_size(dim, interp_order=interp, cells=cells)
-    dl = [1.0 / v for v in cells]
-    b0 = [[3] * dim, [8] * dim]
+
     args = {
         "interp_order": interp,
         "smallest_patch_size": smallest_patch_size,
         "largest_patch_size": [20] * dim,
         "time_step_nbr": 1000,
-        "final_time": 1.0,
+        "time_step": 0.001,
         "boundary_types": ["periodic"] * dim,
         "cells": cells,
-        "dl": dl,
-        "refinement_boxes": {"L0": {"B0": b0}},
+        "dl": [1.0 / v for v in cells],
         "refined_particle_nbr": valid_refined_particle_nbr[dim][interp][0],
         "diag_options": {},
         "nesting_buffer": 0,
-        "strict": True,
+        # "strict": True,
     }
-    for k, v in kwargs.items():
-        if k in args:
-            args[k] = v
+    args.update(deepcopy(kwargs))
+
+    if "refinement" not in kwargs and "refinement_boxes" not in kwargs:
+        b0 = [[3] * dim, [8] * dim]
+        args["refinement_boxes"] = {"L0": {"B0": b0}}
 
+    print("merged", args)
     return args
 
 
@@ -139,20 +145,20 @@ def makeBasicModel(extra_pops={}):
     )
 
 
-def populate_simulation(dim, interp, **input):
+def populate_simulation(ndim=1, interp=1, model_fn=None, diags_fn=None, **simInput):
     ph.global_vars.sim = None
-    simulation = ph.Simulation(**basicSimulatorArgs(dim, interp, **input))
+    simulation = ph.Simulation(**basicSimulatorArgs(ndim, interp, **simInput))
     extra_pops = {}
-    if "populations" in input:
-        for pop, vals in input["populations"].items():
+    if "populations" in simInput:
+        for pop, vals in simInput["populations"].items():
             extra_pops[pop] = defaultPopulationSettings()
             extra_pops[pop].update(vals)
 
-    model = makeBasicModel(extra_pops)
-    if "diags_fn" in input:
-        input["diags_fn"](model)
+    model = model_fn() if model_fn else makeBasicModel(extra_pops)
+    if diags_fn:
+        diags_fn(model)
 
-    ElectronModel(closure="isothermal", Te=0.12)
+    ph.ElectronModel(closure="isothermal", Te=0.12)
 
     return simulation