diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index bcc78a7d7..b869dddca 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -44,6 +44,7 @@ serialize as serialize_sim, deserialize as deserialize_sim, ) +from .load_balancer import LoadBalancer def NO_GUI(): @@ -86,9 +87,7 @@ def __init__(self, fn): self.fn = fn def __call__(self, *xyz): - args = [] - for i, arg in enumerate(xyz): - args.append(np.asarray(arg)) + args = [np.asarray(arg) for arg in xyz] ret = self.fn(*args) if isinstance(ret, list): ret = np.asarray(ret) @@ -128,6 +127,9 @@ def populateDict(): def add_int(path, val): pp.add_int(path, int(val)) + def add_bool(path, val): + pp.add_bool(path, bool(val)) + def add_double(path, val): pp.add_double(path, float(val)) @@ -182,8 +184,6 @@ def add_vector_int(path, val): add_int("simulation/AMR/tag_buffer", simulation.tag_buffer) - add_string("simulation/AMR/loadbalancing", simulation.loadbalancing) - refinement_boxes = simulation.refinement_boxes def as_paths(rb): @@ -223,6 +223,27 @@ def as_paths(rb): add_double("simulation/algo/ohm/resistivity", simulation.resistivity) add_double("simulation/algo/ohm/hyper_resistivity", simulation.hyper_resistivity) + # load balancer block start + lb = simulation.load_balancer or LoadBalancer(_register=False) + base = "simulation/AMR/loadbalancing" + add_bool(f"{base}/active", lb.active) + add_string(f"{base}/mode", lb.mode) + add_double(f"{base}/tolerance", lb.tol) + + # if mode==nppc, imbalance allowed + add_bool(f"{base}/auto", lb.auto) + add_size_t(f"{base}/next_rebalance", lb.next_rebalance) + add_size_t(f"{base}/max_next_rebalance", lb.max_next_rebalance) + add_size_t( + f"{base}/next_rebalance_backoff_multiplier", + lb.next_rebalance_backoff_multiplier, + ) + + # cadence based values + add_size_t(f"{base}/every", lb.every) + add_bool(f"{base}/on_init", lb.on_init) + # load balancer block end + init_model = simulation.model modelDict = init_model.model_dict @@ -246,12 +267,14 @@ def as_paths(rb): addInitFunction(partinit_path + "thermal_velocity_x", fn_wrapper(d["vthx"])) addInitFunction(partinit_path + "thermal_velocity_y", fn_wrapper(d["vthy"])) addInitFunction(partinit_path + "thermal_velocity_z", fn_wrapper(d["vthz"])) - add_int(partinit_path + "nbr_part_per_cell", d["nbrParticlesPerCell"]) add_double(partinit_path + "charge", d["charge"]) add_string(partinit_path + "basis", "cartesian") if "init" in d and "seed" in d["init"]: pp.add_optional_size_t(partinit_path + "init/seed", d["init"]["seed"]) + add_int(partinit_path + "nbr_part_per_cell", d["nbrParticlesPerCell"]) + add_double(partinit_path + "density_cut_off", d["density_cut_off"]) + add_string("simulation/electromag/name", "EM") add_string("simulation/electromag/electric/name", "E") diff --git a/pyphare/pyphare/pharein/load_balancer.py b/pyphare/pyphare/pharein/load_balancer.py new file mode 100644 index 000000000..ee9ebf72e --- /dev/null +++ b/pyphare/pyphare/pharein/load_balancer.py @@ -0,0 +1,50 @@ +# +# + +from dataclasses import dataclass, field +from . import global_vars as gv + + +@dataclass +class LoadBalancer: + # whether or not load balancing is performed + active: bool = field(default_factory=lambda: False) + + # which way load is assessed + mode: str = field(default_factory=lambda: "nppc") + + # acceptable imbalance essentially + tol: float = field(default_factory=lambda: 0.05) + + # whether to rebalance/check imbalance on init + on_init: bool = field(default_factory=lambda: True) + + # if auto, other values are not used if active + auto: bool = field(default_factory=lambda: True) + next_rebalance_backoff_multiplier: int = field(default_factory=lambda: 2) + next_rebalance: int = field(default_factory=lambda: 200) + max_next_rebalance: int = field(default_factory=lambda: 1000) + + # if !auto these values are used if active + every: int = field(default_factory=lambda: 1) + + # internal, allows not registering object for default init + _register: bool = field(default_factory=lambda: True) + + def __post_init__(self): + allowed_modes = [ + "nppc", # count particles per rank + "homogeneous", # count cells per rank + ] + + if self.mode not in allowed_modes: + raise RuntimeError(f"LoadBalancer mode '{self.mode}' is not valid") + + if self._register: + if not gv.sim: + raise RuntimeError( + f"LoadBalancer cannot be registered as no simulation exists" + ) + if gv.sim.load_balancer: + raise RuntimeError(f"LoadBalancer is already registered to simulation") + gv.sim.load_balancer = self diff --git a/pyphare/pyphare/pharein/maxwellian_fluid_model.py b/pyphare/pyphare/pharein/maxwellian_fluid_model.py index ecb7ee7ff..0063bc866 100644 --- a/pyphare/pyphare/pharein/maxwellian_fluid_model.py +++ b/pyphare/pyphare/pharein/maxwellian_fluid_model.py @@ -71,6 +71,7 @@ def add_population( vthy=None, vthz=None, init={}, + density_cut_off=1e-16, ): """ add a particle population to the current model @@ -122,6 +123,7 @@ def add_population( "vthz": vthz, "nbrParticlesPerCell": nbr_part_per_cell, "init": init, + "density_cut_off": density_cut_off, } } diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 448376489..44d85806c 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -121,6 +121,8 @@ def check_time(**kwargs): and "time_step" not in kwargs ) + start_time = kwargs.get("restart_options", {}).get("restart_time", 0) + if final_and_dt: time_step_nbr = int(kwargs["final_time"] / kwargs["time_step"]) time_step = kwargs["final_time"] / time_step_nbr @@ -139,7 +141,11 @@ def check_time(**kwargs): + " or 'final_time' and 'time_step_nbr'" ) - return time_step_nbr, time_step, kwargs.get("final_time", time_step * time_step_nbr) + return ( + time_step_nbr, + time_step, + kwargs.get("final_time", start_time + time_step * time_step_nbr), + ) # ------------------------------------------------------------------------------ @@ -831,6 +837,7 @@ def __init__(self, **kwargs): self.diagnostics = {} self.model = None self.electrons = None + self.load_balancer = None # hard coded in C++ MultiPhysicsIntegrator::getMaxFinerLevelDt self.nSubcycles = 4 @@ -846,9 +853,6 @@ def __init__(self, **kwargs): ] validate_restart_options(self) - def final_time(self): - return self.time_step * self.time_step_nbr - def simulation_domain(self): return [dl * n + ori for dl, n, ori in zip(self.dl, self.cells, self.origin)] diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py index a7231bcac..d3a7fef1a 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -10,6 +10,8 @@ from .particles import Particles from ..core.phare_utilities import listify +from dataclasses import dataclass + class PatchData: """ @@ -86,6 +88,9 @@ def __str__(self): def __repr__(self): return self.__str__() + def __eq__(self, that): + return self.field_name == that.field_name and self.dataset[:] == that.dataset[:] + def select(self, box): """ return view of internal data based on overlap of input box @@ -221,6 +226,9 @@ def __getitem__(self, box): def size(self): return self.dataset.size() + def __eq__(self, that): + return self.name == that.name and self.dataset == that.dataset + class Patch: """ @@ -247,6 +255,9 @@ def __str__(self): def __repr__(self): return self.__str__() + def __getitem__(self, key): + return self.patch_datas[key] + def copy(self): """does not copy patchdatas.datasets (see class PatchData)""" from copy import deepcopy @@ -1755,6 +1766,7 @@ def merge_particles(hierarchy): popname = domain_pdata[0].split("_")[0] pdatas[popname + "_particles"] = pdatas[domain_pdata[0]] del pdatas[domain_pdata[0]] + return hierarchy def h5_filename_from(diagInfo): @@ -1792,3 +1804,58 @@ def getPatch(hier, point): print("error : ", k, v) raise RuntimeError("more than one patch found for point") return patches + + +@dataclass +class EqualityReport: + ok: bool + reason: str + + def __bool__(self): + return self.ok + + +def hierarchy_compare(this, that): + if not isinstance(this, PatchHierarchy) or not isinstance(that, PatchHierarchy): + return EqualityReport(False, "class type mismatch") + + if this.ndim != that.ndim or this.domain_box != that.domain_box: + return EqualityReport(False, "dimensional mismatch") + + if this.time_hier.keys() != that.time_hier.keys(): + return EqualityReport(False, "timesteps mismatch") + + for tidx in this.times(): + patch_levels_ref = this.time_hier[tidx] + patch_levels_cmp = that.time_hier[tidx] + + if patch_levels_ref.keys() != patch_levels_cmp.keys(): + return EqualityReport(False, "levels mismatch") + + for level_idx in patch_levels_cmp.keys(): + patch_level_ref = patch_levels_ref[level_idx] + patch_level_cmp = patch_levels_cmp[level_idx] + + for patch_idx in range(len(patch_level_cmp.patches)): + patch_ref = patch_level_ref.patches[patch_idx] + patch_cmp = patch_level_cmp.patches[patch_idx] + + if patch_ref.patch_datas.keys() != patch_cmp.patch_datas.keys(): + print(list(patch_ref.patch_datas.keys())) + print(list(patch_cmp.patch_datas.keys())) + return EqualityReport(False, "data keys mismatch") + + for patch_data_key in patch_ref.patch_datas.keys(): + patch_data_ref = patch_ref.patch_datas[patch_data_key] + patch_data_cmp = patch_cmp.patch_datas[patch_data_key] + + if patch_data_cmp != patch_data_ref: + return EqualityReport( + False, + "data mismatch: " + + type(patch_data_cmp).__name__ + + " " + + type(patch_data_ref).__name__, + ) + + return EqualityReport(True, "OK") diff --git a/pyphare/pyphare/pharesee/particles.py b/pyphare/pyphare/pharesee/particles.py index fbf1422d2..f831545ed 100644 --- a/pyphare/pyphare/pharesee/particles.py +++ b/pyphare/pyphare/pharesee/particles.py @@ -272,3 +272,32 @@ def _arg_sort(particles): if particles.ndim == 3: z1 = particles.iCells[:, 2] + particles.deltas[:, 2] return np.argsort(np.sqrt((x1**2 + y1**2 + z1**2)) / (x1 / y1 / z1)) + + +def single_patch_per_level_per_pop_from(hier, only_keep_L0=True): # dragons + for tidx in hier.times(): + if only_keep_L0: + hier.time_hier[tidx] = {0: hier.time_hier[tidx][0]} + + patch_levels = hier.time_hier[tidx] + + for level_idx in patch_levels.keys(): + patch_level = patch_levels[level_idx] + + patch0 = patch_level.patches[0] + particles = {} # str:[] + + for key in patch0.patch_datas.keys(): + if isinstance(patch0[key].dataset, Particles): + particles[key] = [] + + for patch in patch_level.patches: + for patch_data_key in patch.patch_datas.keys(): + particles[key] += [patch[patch_data_key].dataset] + + for key in particles.keys(): + patch0[key].dataset = aggregate(particles[key]) + + patch_levels[level_idx].patches = [patch0] # just one patch + + return hier diff --git a/pyphare/pyphare/simulator/simulator.py b/pyphare/pyphare/simulator/simulator.py index 8e7f0110f..50bca812f 100644 --- a/pyphare/pyphare/simulator/simulator.py +++ b/pyphare/pyphare/simulator/simulator.py @@ -1,3 +1,4 @@ +import datetime import atexit import time as timem import numpy as np @@ -29,6 +30,28 @@ def startMPI(): life_cycles["samrai"] = cpp_lib().SamraiLifeCycle() +def print_rank0(*args, **kwargs): + from pyphare.cpp import cpp_lib + + if cpp_lib().mpi_rank() == 0: + print(*args, **kwargs) + + +def plot_timestep_time(timestep_times): + from pyphare.cpp import cpp_lib + + if cpp_lib().mpi_rank() == 0: + import matplotlib.pyplot as plt + + fig, ax = plt.subplots() + ax.plot(timestep_times) + plt.ylabel("timestep time") + plt.xlabel("timestep") + fig.savefig("timestep_times.png") + + cpp_lib().mpi_barrier() + + class Simulator: def __init__(self, simulation, auto_dump=True, **kwargs): import pyphare.pharein as ph @@ -107,8 +130,7 @@ def _throw(self, e): import sys from pyphare.cpp import cpp_lib - if cpp_lib().mpi_rank() == 0: - print(e) + print_rank0(e) sys.exit(1) def advance(self, dt=None): @@ -134,7 +156,7 @@ def times(self): self.timeStep(), ) - def run(self): + def run(self, plot_times=False): from pyphare.cpp import cpp_lib self._check_init() @@ -153,8 +175,11 @@ def run(self): out = f"t = {t:8.5f} - {ticktock:6.5f}sec - total {np.sum(perf):7.4}sec" print(out, end=self.print_eol) - print("mean advance time = {}".format(np.mean(perf))) - print("total advance time = {}".format(np.sum(perf))) + print_rank0(f"mean advance time = {np.mean(perf)}") + print_rank0(f"total advance time = {datetime.timedelta(seconds=np.sum(perf))}") + + if plot_times: + plot_timestep_time(perf) return self.reset() diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_3d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_3d.py new file mode 100644 index 000000000..b4c22a01c --- /dev/null +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_3d.py @@ -0,0 +1 @@ +test_geometry_3d.py diff --git a/res/cmake/def.cmake b/res/cmake/def.cmake index 2d68b4c96..9bf260756 100644 --- a/res/cmake/def.cmake +++ b/res/cmake/def.cmake @@ -142,6 +142,7 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests endif() function(set_exe_paths_ binary) + set_property(TEST ${binary} PROPERTY ENVIRONMENT "PATH=$ENV{PATH}") set_property(TEST ${binary} PROPERTY ENVIRONMENT "PYTHONPATH=${PHARE_PYTHONPATH}") # ASAN detects leaks by default, even in system/third party libraries set_property(TEST ${binary} APPEND PROPERTY ENVIRONMENT "ASAN_OPTIONS=detect_leaks=0") @@ -174,24 +175,24 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests function(add_no_mpi_python3_test name file directory) if(NOT testMPI OR (testMPI AND forceSerialTests)) - add_test(NAME py3_${name} COMMAND python3 -u ${file} WORKING_DIRECTORY ${directory}) + add_test(NAME py3_${name} COMMAND ${Python_EXECUTABLE} -u ${file} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${name}) endif() endfunction(add_no_mpi_python3_test) if(testMPI) function(add_phare_test binary directory) - add_test(NAME ${binary} COMMAND mpirun -n ${PHARE_MPI_PROCS} ./${binary} WORKING_DIRECTORY ${directory}) + add_test(NAME ${binary} COMMAND ${MPIEXEC_EXECUTABLE} -n ${PHARE_MPI_PROCS} ./${binary} WORKING_DIRECTORY ${directory}) add_phare_test_(${binary} ${directory}) endfunction(add_phare_test) function(add_python3_test name file directory) - add_test(NAME py3_${name} COMMAND mpirun -n ${PHARE_MPI_PROCS} python3 -u ${file} WORKING_DIRECTORY ${directory}) + add_test(NAME py3_${name} COMMAND ${MPIEXEC_EXECUTABLE} -n ${PHARE_MPI_PROCS} ${Python_EXECUTABLE} -u ${file} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${name}) endfunction(add_python3_test) function(add_mpi_python3_test N name file directory) - add_test(NAME py3_${name}_mpi_n_${N} COMMAND mpirun -n ${N} python3 ${file} WORKING_DIRECTORY ${directory}) + add_test(NAME py3_${name}_mpi_n_${N} COMMAND ${MPIEXEC_EXECUTABLE} -n ${N} ${Python_EXECUTABLE} ${file} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${name}_mpi_n_${N}) endfunction(add_mpi_python3_test) @@ -248,13 +249,13 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests if(${N} EQUAL 1) add_test( NAME py3_${target} - COMMAND python3 -u ${file} ${CLI_ARGS} + COMMAND ${Python_EXECUTABLE} -u ${file} ${CLI_ARGS} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${target}) else() add_test( NAME py3_${target}_mpi_n_${N} - COMMAND mpirun -n ${N} python3 -u ${file} ${CLI_ARGS} + COMMAND ${MPIEXEC_EXECUTABLE} -n ${N} ${Python_EXECUTABLE} -u ${file} ${CLI_ARGS} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${target}_mpi_n_${N}) endif() @@ -272,7 +273,7 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests function(phare_python3_exec level target file directory) if(${level} GREATER_EQUAL ${PHARE_EXEC_LEVEL_MIN} AND ${level} LESS_EQUAL ${PHARE_EXEC_LEVEL_MAX}) string (REPLACE ";" " " CLI_ARGS "${ARGN}") - add_test(NAME py3_${target} COMMAND python3 -u ${file} ${CLI_ARGS} WORKING_DIRECTORY ${directory}) + add_test(NAME py3_${target} COMMAND ${Python_EXECUTABLE} -u ${file} ${CLI_ARGS} WORKING_DIRECTORY ${directory}) set_exe_paths_(py3_${target}) endif() endfunction(phare_python3_exec) diff --git a/res/cmake/dep/samrai.cmake b/res/cmake/dep/samrai.cmake index 303861da5..beda4895f 100644 --- a/res/cmake/dep/samrai.cmake +++ b/res/cmake/dep/samrai.cmake @@ -3,6 +3,7 @@ find_package(MPI REQUIRED COMPONENTS C) # some linkers (like mold) don't use default library paths get_filename_component(MPI_LIBRARY_PATH ${MPI_LIBRARY} DIRECTORY) + find_package(SAMRAI CONFIG QUIET) if (NOT SAMRAI_FOUND) message("SAMRAI NOT FOUND") diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index 18ecbca9b..187fecbef 100644 --- a/src/amr/level_initializer/hybrid_level_initializer.hpp +++ b/src/amr/level_initializer/hybrid_level_initializer.hpp @@ -52,24 +52,24 @@ namespace solver auto& hybridModel = static_cast(model); auto& level = amr_types::getLevel(*hierarchy, levelNumber); - auto& hybMessenger = dynamic_cast(messenger); + auto& hybMessenger = dynamic_cast(messenger); + bool const isRegriddingL0 = isRegridding and levelNumber == 0; - if (isRootLevel(levelNumber)) + if (isRegridding) { - PHARE_LOG_START(3, "hybridLevelInitializer::initialize : root level init"); - model.initialize(level); - messenger.fillRootGhosts(model, level, initDataTime); - PHARE_LOG_STOP(3, "hybridLevelInitializer::initialize : root level init"); + std::cout << "regriding level " << levelNumber << "\n"; + PHARE_LOG_START(5, "hybridLevelInitializer::initialize : regriding block"); + messenger.regrid(hierarchy, levelNumber, oldLevel, model, initDataTime); + PHARE_LOG_STOP(5, "hybridLevelInitializer::initialize : regriding block"); } - else { - if (isRegridding) + if (isRootLevel(levelNumber)) { - std::cout << "regriding level " << levelNumber << "\n"; - PHARE_LOG_START(3, "hybridLevelInitializer::initialize : regriding block"); - messenger.regrid(hierarchy, levelNumber, oldLevel, model, initDataTime); - PHARE_LOG_STOP(3, "hybridLevelInitializer::initialize : regriding block"); + PHARE_LOG_START(5, "hybridLevelInitializer::initialize : root level init"); + model.initialize(level); + messenger.fillRootGhosts(model, level, initDataTime); + PHARE_LOG_STOP(5, "hybridLevelInitializer::initialize : root level init"); } else { @@ -107,8 +107,9 @@ namespace solver // it seems SAMRAI does not call timeInterpolate() at this point although // both moments and J need time interpolation. It probably knows that // we are at a sync time across levels and that the time interpolation - // is not needed. But is still seems to use the messenger tempoeraries like + // is not needed. But is still seems to use the messenger temporaries like // NiOld etc. so prepareStep() must be called, see end of the function. + // - TODO more better comment(s) hybMessenger.fillIonMomentGhosts(hybridModel.state.ions, level, initDataTime); @@ -116,41 +117,43 @@ namespace solver // via Ampere and Ohm // this only needs to be done for the root level // since otherwise initLevel has done it already - - if (isRootLevel(levelNumber)) - { - auto& B = hybridModel.state.electromag.B; - auto& J = hybridModel.state.J; - - for (auto& patch : level) + // TODO NICO comment! E is regridded, we only needed J for E + if (!isRegriddingL0) + if (isRootLevel(levelNumber)) { - auto _ = hybridModel.resourcesManager->setOnPatch(*patch, B, J); - auto layout = PHARE::amr::layoutFromPatch(*patch); - auto __ = core::SetLayout(&layout, ampere_); - ampere_(B, J); - - hybridModel.resourcesManager->setTime(J, *patch, 0.); + auto& B = hybridModel.state.electromag.B; + auto& J = hybridModel.state.J; + + for (auto& patch : level) + { + auto _ = hybridModel.resourcesManager->setOnPatch(*patch, B, J); + auto layout = PHARE::amr::layoutFromPatch(*patch); + auto __ = core::SetLayout(&layout, ampere_); + ampere_(B, J); + + hybridModel.resourcesManager->setTime(J, *patch, 0.); + } + hybMessenger.fillCurrentGhosts(J, levelNumber, 0.); + + auto& electrons = hybridModel.state.electrons; + auto& E = hybridModel.state.electromag.E; + + for (auto& patch : level) + { + auto layout = PHARE::amr::layoutFromPatch(*patch); + auto _ + = hybridModel.resourcesManager->setOnPatch(*patch, B, E, J, electrons); + electrons.update(layout); + auto& Ve = electrons.velocity(); + auto& Ne = electrons.density(); + auto& Pe = electrons.pressure(); + auto __ = core::SetLayout(&layout, ohm_); + ohm_(Ne, Ve, Pe, B, J, E); + hybridModel.resourcesManager->setTime(E, *patch, 0.); + } + + hybMessenger.fillElectricGhosts(E, levelNumber, 0.); } - hybMessenger.fillCurrentGhosts(J, levelNumber, 0.); - - auto& electrons = hybridModel.state.electrons; - auto& E = hybridModel.state.electromag.E; - - for (auto& patch : level) - { - auto layout = PHARE::amr::layoutFromPatch(*patch); - auto _ = hybridModel.resourcesManager->setOnPatch(*patch, B, E, J, electrons); - electrons.update(layout); - auto& Ve = electrons.velocity(); - auto& Ne = electrons.density(); - auto& Pe = electrons.pressure(); - auto __ = core::SetLayout(&layout, ohm_); - ohm_(Ne, Ve, Pe, B, J, E); - hybridModel.resourcesManager->setTime(E, *patch, 0.); - } - - hybMessenger.fillElectricGhosts(E, levelNumber, 0.); - } // quantities have been computed on the level,like the moments and J // that we later in the code need to get on level ghost nodes via diff --git a/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp b/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp index ab93c7ea2..0db1e7c03 100644 --- a/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp +++ b/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp @@ -3,7 +3,7 @@ #ifndef PHARE_CONCRERTE_LOAD_BALANCER_HYBRID_STRATEGY_HOMOGENEOUS_HPP #define PHARE_CONCRERTE_LOAD_BALANCER_HYBRID_STRATEGY_HOMOGENEOUS_HPP -#include + #include #include @@ -15,7 +15,6 @@ - namespace PHARE::amr { template @@ -24,14 +23,16 @@ class ConcreteLoadBalancerHybridStrategyHomogeneous : public LoadBalancerHybridS public: using HybridModel = typename PHARE_T::HybridModel_t; using gridlayout_type = typename HybridModel::gridlayout_type; + using amr_types = typename HybridModel::amr_types; + using level_t = typename amr_types::level_t; + using cell_data_t = SAMRAI::pdat::CellData; ConcreteLoadBalancerHybridStrategyHomogeneous(int const id) : id_{id} { } - void compute(SAMRAI::hier::PatchLevel& level, - PHARE::solver::IPhysicalModel& model) override; + void compute(level_t& level, PHARE::solver::IPhysicalModel& model) override; private: @@ -42,30 +43,14 @@ class ConcreteLoadBalancerHybridStrategyHomogeneous : public LoadBalancerHybridS template void ConcreteLoadBalancerHybridStrategyHomogeneous::compute( - SAMRAI::hier::PatchLevel& level, PHARE::solver::IPhysicalModel& model) + level_t& level, PHARE::solver::IPhysicalModel& model) { - static auto constexpr dimension = HybridModel::dimension; - auto& hybridModel = dynamic_cast(model); - auto& resourcesManager = hybridModel.resourcesManager; - auto& ions = hybridModel.state.ions; - - bool constexpr c_ordering = false; - for (auto& patch : level) { - auto const& layout = layoutFromPatch(*patch); - - auto patch_data_lb - = dynamic_cast*>(patch->getPatchData(this->id_).get()); + auto const& layout = layoutFromPatch(*patch); + auto patch_data_lb = dynamic_cast(patch->getPatchData(this->id_).get()); auto load_balancer_val = patch_data_lb->getPointer(); - const auto& box = patch->getBox(); - - auto lb_view = core::NdArrayView(load_balancer_val, - layout.nbrCells()); - - auto _ = resourcesManager->setOnPatch(*patch, ions); - // The view on "load_balancer_val" do not have any ghost cells, meaning that this patch // data lies on a box defind by the nbr of cells in the physical domain // As the nbr of ghost cells in the box associated to the load balancer patch data is @@ -73,64 +58,10 @@ void ConcreteLoadBalancerHybridStrategyHomogeneous::compute( // to get the AMRLocatStart and AMRLocalEnd. // Then, we build "by hand" the local index for the "lb_view" considering that the nbr // of ghost cell for this patch data is null. - auto amr_box = layout.AMRBox(); // The lb_view is a CellData, meaning that it is dual as the index of an amr box - auto amr_start = amr_box.lower; - auto amr_end = amr_box.upper; - - // nbr of cells in the physical domain - auto nbrCells = layout.nbrCells(); - - - - - if constexpr (dimension == 1) - { - for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - auto amr_point{core::Point{i_amr}}; - auto amr_index = amr_point.template toArray(); - - lb_view(i_loc) = 1; - } - } - - - - if constexpr (dimension == 2) - { - for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - for (auto j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - { - auto amr_point{core::Point{i_amr, j_amr}}; - auto amr_index = amr_point.template toArray(); - - lb_view(i_loc, j_loc) = 1; - std::cout << lb_view(i_loc, j_loc) << std::endl; - } - } - } - - - if constexpr (dimension == 3) - { - for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - for (auto j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - { - for (auto k_loc = 0, k_amr = amr_start[2]; k_loc < nbrCells[2]; - ++k_loc, ++k_amr) - { - auto amr_point{core::Point{i_amr, j_amr, k_amr}}; - auto amr_index = amr_point.template toArray(); - lb_view(i_loc, j_loc, k_loc) = 1; - } - } - } - } + std::fill(load_balancer_val, load_balancer_val + core::product(layout.nbrCells()), 1); } } diff --git a/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_nppc.hpp b/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_nppc.hpp index 6965a0f1a..d97a77190 100644 --- a/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_nppc.hpp +++ b/src/amr/load_balancing/concrete_load_balancer_hybrid_strategy_nppc.hpp @@ -8,11 +8,14 @@ #include +#include "core/logger.hpp" +#include "core/utilities/types.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" + +#include "amr/types/amr_types.hpp" #include "amr/load_balancing/load_balancer_hybrid_strategy.hpp" #include "amr/physical_models/physical_model.hpp" -#include "amr/types/amr_types.hpp" #include "amr/resources_manager/amr_utils.hpp" -#include "core/data/ndarray/ndarray_vector.hpp" @@ -25,14 +28,16 @@ class ConcreteLoadBalancerHybridStrategyNPPC : public LoadBalancerHybridStrategy public: using HybridModel = typename PHARE_T::HybridModel_t; using gridlayout_type = typename HybridModel::gridlayout_type; + using amr_types = typename HybridModel::amr_types; + using level_t = typename amr_types::level_t; + using cell_data_t = SAMRAI::pdat::CellData; ConcreteLoadBalancerHybridStrategyNPPC(int const id) : id_{id} { } - void compute(SAMRAI::hier::PatchLevel& level, - PHARE::solver::IPhysicalModel& model) override; + void compute(level_t& level, PHARE::solver::IPhysicalModel& model) override; private: @@ -43,29 +48,22 @@ class ConcreteLoadBalancerHybridStrategyNPPC : public LoadBalancerHybridStrategy template void ConcreteLoadBalancerHybridStrategyNPPC::compute( - SAMRAI::hier::PatchLevel& level, PHARE::solver::IPhysicalModel& model) + level_t& level, PHARE::solver::IPhysicalModel& model) { - static auto constexpr dimension = HybridModel::dimension; - auto& hybridModel = dynamic_cast(model); - auto& resourcesManager = hybridModel.resourcesManager; - auto& ions = hybridModel.state.ions; + bool static constexpr c_ordering = false; + auto static constexpr dimension = HybridModel::dimension; - bool constexpr c_ordering = false; + auto& hybridModel = dynamic_cast(model); + auto& resourcesManager = hybridModel.resourcesManager; + auto& ions = hybridModel.state.ions; for (auto& patch : level) { - auto const& layout = layoutFromPatch(*patch); - - auto patch_data_lb - = dynamic_cast*>(patch->getPatchData(this->id_).get()); + auto const& layout = layoutFromPatch(*patch); + auto patch_data_lb = dynamic_cast(patch->getPatchData(this->id_).get()); auto load_balancer_val = patch_data_lb->getPointer(); - - const auto& box = patch->getBox(); - - auto lb_view = core::NdArrayView(load_balancer_val, - layout.nbrCells()); - - auto _ = resourcesManager->setOnPatch(*patch, ions); + auto lb_view = core::make_array_view(load_balancer_val, layout.nbrCells()); + auto _ = resourcesManager->setOnPatch(*patch, ions); // The view on "load_balancer_val" do not have any ghost cells, meaning that this patch // data lies on a box defind by the nbr of cells in the physical domain @@ -74,98 +72,27 @@ void ConcreteLoadBalancerHybridStrategyNPPC::compute( // to get the AMRLocatStart and AMRLocalEnd. // Then, we build "by hand" the local index for the "lb_view" considering that the nbr // of ghost cell for this patch data is null. - auto amr_box = layout.AMRBox(); // The lb_view is a CellData, meaning that it is dual as the index of an amr box - auto amr_start = amr_box.lower; - auto amr_end = amr_box.upper; - - // nbr of cells in the physical domain - auto nbrCells = layout.nbrCells(); - - - - - if constexpr (dimension == 1) - { - for (std::uint32_t i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - auto amr_point{core::Point{i_amr}}; - auto amr_index = amr_point.template toArray(); - - std::size_t nbr = 0; - - for (auto& pop : ions) - { - const auto& domainParticles = pop.domainParticles(); - - nbr += domainParticles.nbr_particles_in(amr_index); - } - - lb_view(i_loc) = nbr; - } - } - - - - if constexpr (dimension == 2) - { - for (std::uint32_t i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - for (std::uint32_t j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - { - auto amr_point{core::Point{i_amr, j_amr}}; - auto amr_index = amr_point.template toArray(); - - std::size_t nbr = 0; - - for (auto& pop : ions) - { - const auto& domainParticles = pop.domainParticles(); - - nbr += domainParticles.nbr_particles_in(amr_index); - } - - lb_view(i_loc, j_loc) = nbr; - } - } - } - - - if constexpr (dimension == 3) - { - for (std::uint32_t i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - { - for (std::uint32_t j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - { - for (std::uint32_t k_loc = 0, k_amr = amr_start[2]; k_loc < nbrCells[2]; - ++k_loc, ++k_amr) - { - auto amr_point{core::Point{i_amr, j_amr, k_amr}}; - auto amr_index = amr_point.template toArray(); - - std::size_t nbr = 0; - - for (auto& pop : ions) - { - const auto& domainParticles = pop.domainParticles(); - - nbr += domainParticles.nbr_particles_in(amr_index); - } - - lb_view(i_loc, j_loc, k_loc) = nbr; - } - } - } - } + core::Box local_box{ + core::Point{core::ConstArray()}, + core::Point{ + core::generate([](auto const& nCell) { return nCell - 1; }, layout.nbrCells())}}; + + auto amr_iter = layout.AMRBox().begin(); + auto lcl_iter = local_box.begin(); + + for (; lcl_iter != local_box.end(); ++amr_iter, ++lcl_iter) + lb_view(*lcl_iter) = core::sum_from(ions, [&](auto const& pop) { + return pop.domainParticles().nbr_particles_in((*amr_iter).toArray()); + }); } // TODO here, we have the lb_view value correctly set on all patches. we also know the id_ - // so this is where we should call setWorkloadPatchDataIndex... which is a method of the CascadePartitioner - // lb_view is a local container containing the datz - // the loadbalancezrmanager knows the id, as well as the loadbalancerestimator - // and the loadbalancerestimator is a cascadpartotioner - + // so this is where we should call setWorkloadPatchDataIndex... which is a method of the + // CascadePartitioner lb_view is a local container containing the data the LoadBalancerNanager + // knows the id, as well as the LoadBalancerEstimator and the LoadBalancerEstimator is a + // CascadePartitioner } } // namespace PHARE::amr diff --git a/src/amr/load_balancing/load_balancer_details.hpp b/src/amr/load_balancing/load_balancer_details.hpp new file mode 100644 index 000000000..5dcd85911 --- /dev/null +++ b/src/amr/load_balancing/load_balancer_details.hpp @@ -0,0 +1,57 @@ +#ifndef PHARE_AMR_LOAD_BALANCER_LOAD_BALANCER_DETAILS_HPP +#define PHARE_AMR_LOAD_BALANCER_LOAD_BALANCER_DETAILS_HPP + +#include +#include + +#include "core/logger.hpp" +#include "initializer/data_provider.hpp" + +namespace PHARE::amr +{ +struct LoadBalancerDetailsDefaults +{ + // is instantiated in constexpr context so following fields are constexpr in that context + std::size_t next_rebalance_backoff_multiplier = 2; + std::size_t next_rebalance = 200; + std::size_t max_next_rebalance = 1000; + double tolerance = .05; +}; + +struct LoadBalancerDetails +{ + LoadBalancerDetailsDefaults constexpr static defaults{}; + + bool const active = false; + bool const automatic = false; + bool const on_init = false; + + std::size_t const every = 0; + std::string const mode; + + double const tolerance = defaults.tolerance; + + std::size_t next_rebalance_backoff_multiplier = defaults.next_rebalance_backoff_multiplier; + std::size_t next_rebalance = defaults.next_rebalance; + std::size_t max_next_rebalance = defaults.max_next_rebalance; + + LoadBalancerDetails static FROM(initializer::PHAREDict const& dict) + { + return { + cppdict::get_value(dict, "active", false), + cppdict::get_value(dict, "auto", false), + cppdict::get_value(dict, "on_init", false), + cppdict::get_value(dict, "every", std::size_t{0}), + cppdict::get_value(dict, "mode", std::string{"nppc"}), + cppdict::get_value(dict, "tol", defaults.tolerance), + cppdict::get_value(dict, "next_rebalance_backoff_multiplier", + defaults.next_rebalance_backoff_multiplier), + cppdict::get_value(dict, "next_rebalance", defaults.next_rebalance), + cppdict::get_value(dict, "max_next_rebalance", defaults.max_next_rebalance), + }; + } +}; + +} // namespace PHARE::amr + +#endif /* PHARE_AMR_LOAD_BALANCER_LOAD_BALANCER_DETAILS_HPP */ diff --git a/src/amr/load_balancing/load_balancer_estimator.hpp b/src/amr/load_balancing/load_balancer_estimator.hpp index 8b0dc8ef3..780edee1a 100644 --- a/src/amr/load_balancing/load_balancer_estimator.hpp +++ b/src/amr/load_balancing/load_balancer_estimator.hpp @@ -25,7 +25,7 @@ class LoadBalancerEstimator virtual ~LoadBalancerEstimator() = default; virtual void estimate(SAMRAI::hier::PatchLevel& level, - PHARE::solver::IPhysicalModel& model) + solver::IPhysicalModel& model) = 0; diff --git a/src/amr/load_balancing/load_balancer_estimator_hybrid.hpp b/src/amr/load_balancing/load_balancer_estimator_hybrid.hpp index dc87cfa9b..fd551e9a3 100644 --- a/src/amr/load_balancing/load_balancer_estimator_hybrid.hpp +++ b/src/amr/load_balancing/load_balancer_estimator_hybrid.hpp @@ -3,172 +3,44 @@ #include #include -#include -#include "load_balancer_estimator.hpp" -// #include "amr/resources_manager/amr_utils.hpp" + #include "amr/physical_models/physical_model.hpp" -// #include "core/data/ndarray/ndarray_vector.hpp" -// #include "amr/resources_manager/amr_utils.hpp" -// #include "core/utilities/point/point.hpp" + +#include "load_balancer_estimator.hpp" #include "load_balancer_hybrid_strategy.hpp" #include "load_balancer_hybrid_strategy_factory.hpp" - - namespace PHARE::amr { template class LoadBalancerEstimatorHybrid : public LoadBalancerEstimator { - using HybridModel = typename PHARE_T::HybridModel_t; - using gridlayout_type = typename HybridModel::gridlayout_type; - - + using HybridModel = typename PHARE_T::HybridModel_t; + using amr_types = typename HybridModel::amr_types; + using level_t = typename amr_types::level_t; public: - // LoadBalancerEstimatorHybrid(int const id) LoadBalancerEstimatorHybrid(std::string strategy_name, int const id) : LoadBalancerEstimator{id} , strat_{LoadBalancerHybridStrategyFactory::create(strategy_name, id)} { } - ~LoadBalancerEstimatorHybrid() = default; // the implementation of a virtual class NEEDS a dtor + // the implementation of a virtual class NEEDS a dtor + ~LoadBalancerEstimatorHybrid() = default; - void estimate(SAMRAI::hier::PatchLevel& level, - solver::IPhysicalModel& model) override; + void estimate(level_t& level, solver::IPhysicalModel& model) override + { + strat_->compute(level, model); + } private: std::unique_ptr> strat_; }; - - -template -inline void -LoadBalancerEstimatorHybrid::estimate(SAMRAI::hier::PatchLevel& level, - solver::IPhysicalModel& model) -{ - strat_->compute(level, model); - - // static auto constexpr dimension = HybridModel::dimension; - // auto& hybridModel = dynamic_cast(model); - // auto& resourcesManager = hybridModel.resourcesManager; - // auto& ions = hybridModel.state.ions; - - // bool constexpr c_ordering = false; - - // for (auto& patch : level) - // { - // auto const& layout = layoutFromPatch(*patch); - - // auto patch_data_lb - // = - // dynamic_cast*>(patch->getPatchData(this->id_).get()); - // auto load_balancer_val = patch_data_lb->getPointer(); - - // const auto& box = patch->getBox(); - - // auto lb_view = core::NdArrayView(load_balancer_val, - // layout.nbrCells()); - - // auto _ = resourcesManager->setOnPatch(*patch, ions); - - // // The view on "load_balancer_val" do not have any ghost cells, meaning that this patch - // // data lies on a box defind by the nbr of cells in the physical domain - // // As the nbr of ghost cells in the box associated to the load balancer patch data is - // // null, we hence loop on an AMR index (then needing to firstly get the AMR box) - // // to get the AMRLocatStart and AMRLocalEnd. - // // Then, we build "by hand" the local index for the "lb_view" considering that the nbr - // // of ghost cell for this patch data is null. - // auto amr_box = layout.AMRBox(); - - // // The lb_view is a CellData, meaning that it is dual as the index of an amr box - // auto amr_start = amr_box.lower; - // auto amr_end = amr_box.upper; - - // // nbr of cells in the physical domain - // auto nbrCells = layout.nbrCells(); - - - // if constexpr (dimension == 1) - // { - // for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - // { - // auto amr_point{core::Point{i_amr}}; - // auto amr_index = amr_point.template toArray(); - - // auto nbr = 0; - - // for (auto& pop : ions) - // { - // const auto& domainParticles = pop.domainParticles(); - - // nbr += domainParticles.nbr_particles_in(amr_index); - // } - - // lb_view(i_loc) = nbr; - // } - // } - - - // if constexpr (dimension == 2) - // { - // for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - // { - // for (auto j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - // { - // auto amr_point{core::Point{i_amr, j_amr}}; - // auto amr_index = amr_point.template toArray(); - - // auto nbr = 0; - - // for (auto& pop : ions) - // { - // const auto& domainParticles = pop.domainParticles(); - - // nbr += domainParticles.nbr_particles_in(amr_index); - // } - - // lb_view(i_loc, j_loc) = nbr; - // } - // } - // } - - - // if constexpr (dimension == 3) - // { - // for (auto i_loc = 0, i_amr = amr_start[0]; i_loc < nbrCells[0]; ++i_loc, ++i_amr) - // { - // for (auto j_loc = 0, j_amr = amr_start[1]; j_loc < nbrCells[1]; ++j_loc, ++j_amr) - // { - // for (auto k_loc = 0, k_amr = amr_start[2]; k_loc < nbrCells[2]; - // ++k_loc, ++k_amr) - // { - // auto amr_point{core::Point{i_amr, j_amr, k_amr}}; - // auto amr_index = amr_point.template toArray(); - - // auto nbr = 0; - - // for (auto& pop : ions) - // { - // const auto& domainParticles = pop.domainParticles(); - - // nbr += domainParticles.nbr_particles_in(amr_index); - // } - - // lb_view(i_loc, j_loc, k_loc) = nbr; - // } - // } - // } - // } - // } -} - } // namespace PHARE::amr #endif diff --git a/src/amr/load_balancing/load_balancer_hybrid_strategy.hpp b/src/amr/load_balancing/load_balancer_hybrid_strategy.hpp index 59edef854..a497e4567 100644 --- a/src/amr/load_balancing/load_balancer_hybrid_strategy.hpp +++ b/src/amr/load_balancing/load_balancer_hybrid_strategy.hpp @@ -14,10 +14,14 @@ namespace PHARE::amr template class LoadBalancerHybridStrategy { + using HybridModel = typename PHARE_T::HybridModel_t; + using amr_types = typename HybridModel::amr_types; + using level_t = typename amr_types::level_t; + public: - virtual void compute(SAMRAI::hier::PatchLevel& level, - PHARE::solver::IPhysicalModel& model) - = 0; + virtual ~LoadBalancerHybridStrategy() {} + + virtual void compute(level_t& level, PHARE::solver::IPhysicalModel& model) = 0; }; diff --git a/src/amr/load_balancing/load_balancer_hybrid_strategy_factory.hpp b/src/amr/load_balancing/load_balancer_hybrid_strategy_factory.hpp index f85a74f09..d03f283a9 100644 --- a/src/amr/load_balancing/load_balancer_hybrid_strategy_factory.hpp +++ b/src/amr/load_balancing/load_balancer_hybrid_strategy_factory.hpp @@ -1,4 +1,3 @@ - #ifndef LOAD_BALANCER_HYBRID_STRATEGY_FACTORY_HPP #define LOAD_BALANCER_HYBRID_STRATEGY_FACTORY_HPP @@ -7,7 +6,7 @@ #include "amr/load_balancing/load_balancer_hybrid_strategy.hpp" #include "amr/load_balancing/concrete_load_balancer_hybrid_strategy_nppc.hpp" - +#include "amr/load_balancing/concrete_load_balancer_hybrid_strategy_homogeneous.hpp" namespace PHARE::amr @@ -26,13 +25,15 @@ class LoadBalancerHybridStrategyFactory else if (strat_name == "homogeneous") { - return std::make_unique>(id); + return std::make_unique>(id); } return nullptr; } }; + } // namespace PHARE::amr + #endif diff --git a/src/amr/load_balancing/load_balancer_manager.hpp b/src/amr/load_balancing/load_balancer_manager.hpp index d3accc7bc..c9a021789 100644 --- a/src/amr/load_balancing/load_balancer_manager.hpp +++ b/src/amr/load_balancing/load_balancer_manager.hpp @@ -7,21 +7,19 @@ #include #include -// #include "phare_core.hpp" + #include "initializer/data_provider.hpp" #include "load_balancer_estimator.hpp" -// #include "amr/resources_manager/amr_utils.hpp" -// #include "amr/solvers/solver.hpp" - namespace PHARE::amr { + + template class LoadBalancerManager { public: - // LoadBalancerManager(int const maxLevelNumber) LoadBalancerManager(PHARE::initializer::PHAREDict const& dict) : dim_{SAMRAI::tbox::Dimension{dim}} , loadBalancerVar_{std::make_shared>( @@ -31,19 +29,18 @@ class LoadBalancerManager , id_{variableDatabase_->registerVariableAndContext(loadBalancerVar_, context_, SAMRAI::hier::IntVector::getZero(dim_))} , maxLevelNumber_{dict["simulation"]["AMR"]["max_nbr_levels"].template to()} - // , loadBalancerEstimators_(maxLevelNumber, nullptr){}; - , loadBalancerEstimators_(maxLevelNumber_, nullptr){}; + , loadBalancerEstimators_(maxLevelNumber_){}; ~LoadBalancerManager() { variableDatabase_->removeVariable("LoadBalancerVariable"); }; - int getId() const; + int getId() const { return id_; } void addLoadBalancerEstimator(int const iLevel_min, int const iLevel_max, std::shared_ptr lbe); - void addLoadBalancer(std::unique_ptr loadBalancer) + void setLoadBalancer(std::shared_ptr loadBalancer) { - loadBalancer_ = std::move(loadBalancer); + loadBalancer_ = loadBalancer; loadBalancer_->setWorkloadPatchDataIndex(id_); } @@ -61,27 +58,17 @@ class LoadBalancerManager int const id_; int const maxLevelNumber_; std::vector> loadBalancerEstimators_; - std::unique_ptr loadBalancer_; + std::shared_ptr loadBalancer_; }; template -inline int LoadBalancerManager::getId() const -{ - return id_; -} - - - - -template -inline void -LoadBalancerManager::addLoadBalancerEstimator(int const iLevel_min, int const iLevel_max, - std::shared_ptr lbe) +void LoadBalancerManager::addLoadBalancerEstimator( + int const iLevel_min, int const iLevel_max, std::shared_ptr lbe) { - for (auto ilevel = iLevel_min; ilevel <= iLevel_max; ilevel++) + for (auto ilevel = iLevel_min; ilevel <= iLevel_max; ++ilevel) { loadBalancerEstimators_[ilevel] = lbe; } @@ -91,8 +78,7 @@ LoadBalancerManager::addLoadBalancerEstimator(int const iLevel_min, int con template -inline void LoadBalancerManager::allocate(SAMRAI::hier::Patch& patch, - double const allocateTime) +void LoadBalancerManager::allocate(SAMRAI::hier::Patch& patch, double const allocateTime) { patch.allocatePatchData(id_, allocateTime); } @@ -100,16 +86,15 @@ inline void LoadBalancerManager::allocate(SAMRAI::hier::Patch& patch, template -inline void -LoadBalancerManager::estimate(SAMRAI::hier::PatchLevel& level, - PHARE::solver::IPhysicalModel& model) +void LoadBalancerManager::estimate( + SAMRAI::hier::PatchLevel& level, PHARE::solver::IPhysicalModel& model) { - auto iLevel = level.getLevelNumber(); - auto lbe = loadBalancerEstimators_[iLevel]; - - lbe->estimate(level, model); + if (auto lbe = loadBalancerEstimators_[level.getLevelNumber()]) + lbe->estimate(level, model); } + } // namespace PHARE::amr + #endif diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp index 3bee625e6..cb0d23a26 100644 --- a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp +++ b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp @@ -1,6 +1,7 @@ #ifndef PHARE_HYBRID_HYBRID_MESSENGER_STRATEGY_HPP #define PHARE_HYBRID_HYBRID_MESSENGER_STRATEGY_HPP +#include "core/logger.hpp" #include "core/def/phare_mpi.hpp" #include "SAMRAI/hier/CoarseFineBoundary.h" @@ -109,7 +110,7 @@ namespace amr using DefaultCoarsenOp = BaseCoarsenOp>; public: - static const std::string stratName; + static const inline std::string stratName = "HybridModel-HybridModel"; static constexpr std::size_t rootLevelNumber = 0; @@ -190,6 +191,7 @@ namespace amr patchGhostPartRefiners_.registerLevel(hierarchy, level); + // root level is not initialized with a schedule using coarser level data // so we don't create these schedules if root level // TODO this 'if' may not be OK if L0 is regrided @@ -223,11 +225,15 @@ namespace amr { auto& hybridModel = dynamic_cast(model); auto level = hierarchy->getPatchLevel(levelNumber); + + bool isRegriddingL0 = levelNumber == 0 and oldLevel; + magneticInitRefiners_.regrid(hierarchy, levelNumber, oldLevel, initDataTime); electricInitRefiners_.regrid(hierarchy, levelNumber, oldLevel, initDataTime); domainParticlesRefiners_.regrid(hierarchy, levelNumber, oldLevel, initDataTime); patchGhostPartRefiners_.fill(levelNumber, initDataTime); + // regriding will fill the new level wherever it has points that overlap // old level. This will include its level border points. // These new level border points will thus take values that where previous @@ -236,14 +242,18 @@ namespace amr // Specifically, we need all fine faces to have equal magnetic field and also // equal to that of the shared coarse face. // This means that we now need to fill ghosts and border included - auto& B = hybridModel.state.electromag.B; - auto& E = hybridModel.state.electromag.E; - // magSharedNodesRefiners_.fill(B, levelNumber, initDataTime); - magGhostsRefiners_.fill(B, levelNumber, initDataTime); - // elecSharedNodesRefiners_.fill(E, levelNumber, initDataTime); - elecGhostsRefiners_.fill(E, levelNumber, initDataTime); - fix_magnetic_divergence_(*hierarchy, levelNumber, B); + if (!isRegriddingL0) + { + auto& B = hybridModel.state.electromag.B; + auto& E = hybridModel.state.electromag.E; + // magSharedNodesRefiners_.fill(B, levelNumber, initDataTime); + magGhostsRefiners_.fill(B, levelNumber, initDataTime); + // elecSharedNodesRefiners_.fill(E, levelNumber, initDataTime); + elecGhostsRefiners_.fill(E, levelNumber, initDataTime); + + fix_magnetic_divergence_(*hierarchy, levelNumber, B); + } // we now call only levelGhostParticlesOld.fill() and not .regrid() // regrid() would refine from next coarser in regions of level not overlaping @@ -253,8 +263,14 @@ namespace amr // https://github.com/PHAREHUB/PHARE/issues/604 calling .fill() ensures that // levelGhostParticlesOld particles are filled exclusively from spliting next // coarser domain ones like when a new finest level is created. - lvlGhostPartOldRefiners_.fill(levelNumber, initDataTime); - copyLevelGhostOldToPushable_(*level, model); + + + if (levelNumber != rootLevelNumber) + { + lvlGhostPartOldRefiners_.fill(levelNumber, initDataTime); + copyLevelGhostOldToPushable_(*level, model); + } + // computeIonMoments_(*level, model); // levelGhostNew will be refined in next firstStep @@ -350,7 +366,7 @@ namespace amr void fillIonGhostParticles(IonsT& ions, SAMRAI::hier::PatchLevel& level, double const fillTime) override { - PHARE_LOG_SCOPE(3, "HybridHybridMessengerStrategy::fillIonGhostParticles"); + PHARE_LOG_SCOPE(1, "HybridHybridMessengerStrategy::fillIonGhostParticles"); for (auto patch : level) { @@ -367,7 +383,7 @@ namespace amr /** - * @brief fillIonMomentGhosts works on moment ghost nodes + * @brief fillIonPopMomentGhosts works on moment ghost nodes * * patch border node moments are completed by the deposition of patch ghost * particles for all populations level border nodes are completed by the deposition @@ -377,7 +393,7 @@ namespace amr void fillIonPopMomentGhosts(IonsT& ions, SAMRAI::hier::PatchLevel& level, double const afterPushTime) override { - PHARE_LOG_SCOPE(3, "HybridHybridMessengerStrategy::fillIonMomentGhosts"); + PHARE_LOG_SCOPE(1, "HybridHybridMessengerStrategy::fillIonMomentGhosts"); auto alpha = timeInterpCoef_(afterPushTime, level.getLevelNumber()); if (level.getLevelNumber() > 0 and (alpha < 0 or alpha > 1)) @@ -1086,9 +1102,6 @@ namespace amr CoarsenOperator_ptr magneticCoarseningOp_{std::make_shared()}; }; - template - const std::string HybridHybridMessengerStrategy::stratName - = "HybridModel-HybridModel"; } // namespace amr diff --git a/src/amr/multiphysics_integrator.hpp b/src/amr/multiphysics_integrator.hpp index 6794b0d13..8cd7be130 100644 --- a/src/amr/multiphysics_integrator.hpp +++ b/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" @@ -48,7 +49,6 @@ namespace solver int solverIndex = NOT_SET; int resourcesManagerIndex = NOT_SET; int taggerIndex = NOT_SET; - // int loadBalancerIndex = NOT_SET; std::string messengerName; }; @@ -300,7 +300,7 @@ namespace solver */ void initializeLevelData(std::shared_ptr const& hierarchy, int const levelNumber, double const initDataTime, - bool const /*canBeRefined*/, bool const /*initialTime*/, + bool const canBeRefined, bool const initialTime, std::shared_ptr const& oldLevel = std::shared_ptr(), bool const allocateData = true) override @@ -311,14 +311,15 @@ namespace solver auto& levelInitializer = getLevelInitializer(model.name()); bool const isRegridding = oldLevel != nullptr; - auto level = hierarchy->getPatchLevel(levelNumber); + std::cout << "init level " << levelNumber << " with regriding = " << isRegridding << "\n"; PHARE_LOG_START(3, "initializeLevelData::allocate block"); + if (allocateData) { - for (auto patch : *level) + for (auto patch : *hierarchy->getPatchLevel(levelNumber)) { model.allocate(*patch, initDataTime); solver.allocate(model, *patch, initDataTime); @@ -328,13 +329,20 @@ namespace solver } PHARE_LOG_STOP(3, "initializeLevelData::allocate block"); - if (isRegridding) + + bool const isRegriddingL0 = isRegridding and levelNumber == 0; + + if (isRegriddingL0) + { + messenger.registerLevel(hierarchy, levelNumber); + } + else if (isRegridding) { // regriding the current level has broken schedules for which // this level is the source or destination // we therefore need to rebuild them - auto finestLvlNbr = hierarchy->getFinestLevelNumber(); - auto nextFiner = (levelNumber == finestLvlNbr) ? levelNumber : levelNumber + 1; + auto const finestLvlNbr = hierarchy->getFinestLevelNumber(); + auto nextFiner = (levelNumber == finestLvlNbr) ? levelNumber : levelNumber + 1; for (auto ilvl = levelNumber; ilvl <= nextFiner; ++ilvl) { @@ -350,6 +358,16 @@ namespace solver levelInitializer.initialize(hierarchy, levelNumber, oldLevel, model, messenger, initDataTime, isRegridding); + + if (isRegriddingL0) + { + for (auto ilvl = 1; ilvl <= hierarchy->getFinestLevelNumber(); ++ilvl) + messenger.registerLevel(hierarchy, ilvl); + + solver.onRegrid(); + } + else + load_balancer_manager_->estimate(*hierarchy->getPatchLevel(levelNumber), model); } @@ -374,6 +392,9 @@ namespace solver auto time = dict_["restarts"]["restart_time"].template to(); load_balancer_manager_->allocate(*patch, time); } + // if load balance on restart advance + load_balancer_manager_->estimate(*hierarchy->getPatchLevel(ilvl), + getModel_(ilvl)); } restartInitialized_ = true; } @@ -478,10 +499,12 @@ namespace solver if (regridAdvance) throw std::runtime_error("Error - regridAdvance must be False and is True"); - auto iLevel = level->getLevelNumber(); - std::cout << "advanceLevel " << iLevel << " with dt = " << newTime - currentTime - << " from t = " << currentTime << "to t = " << newTime << "\n"; + + PHARE_LOG_LINE_STR("advanceLevel " << iLevel << " with dt = " << newTime - currentTime + << " from t = " << currentTime + << " to t = " << newTime); + auto& solver = getSolver_(iLevel); auto& model = getModel_(iLevel); auto& fromCoarser = getMessengerWithCoarser_(iLevel); @@ -619,20 +642,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; @@ -887,7 +896,6 @@ namespace solver { auto& descriptor = levelDescriptors_[iLevel]; auto messenger = messengers_[descriptor.messengerName].get(); - auto s = messenger->name(); if (messenger) return *messenger; diff --git a/src/amr/wrappers/integrator.hpp b/src/amr/wrappers/integrator.hpp index 25455e70d..ed03c1769 100644 --- a/src/amr/wrappers/integrator.hpp +++ b/src/amr/wrappers/integrator.hpp @@ -1,8 +1,10 @@ #ifndef INTEGRATOR_HPP #define INTEGRATOR_HPP +#include "core/logger.hpp" #include "core/def/phare_mpi.hpp" +#include #include #include #include @@ -13,7 +15,6 @@ #include #include #include -// #include #include #include #include @@ -21,31 +22,115 @@ #include - #include "initializer/data_provider.hpp" +#include "amr/load_balancing/load_balancer_details.hpp" + namespace PHARE::amr { + template class Integrator { + bool static _is_tagging_refinement(initializer::PHAREDict const& dict) + { + return cppdict::get_value(dict, "simulation/AMR/refinement/tagging/method", + std::string{"none"}) + == std::string{"auto"}; + } + + public: static constexpr std::size_t dimension = _dimension; - double advance(double dt) { return timeRefIntegrator_->advanceHierarchy(dt); } + double advance(double dt) + { + bool rebalance_coarsest_now = _should_rebalance_now(); - void initialize() { timeRefIntegrator_->initializeHierarchy(); } + auto new_time = timeRefIntegrator_->advanceHierarchy(dt, rebalance_coarsest_now); + ++time_step_idx; + return new_time; + } + void initialize() { timeRefIntegrator_->initializeHierarchy(); } Integrator(PHARE::initializer::PHAREDict const& dict, std::shared_ptr hierarchy, std::shared_ptr timeRefLevelStrategy, std::shared_ptr tagAndInitStrategy, - double startTime, double endTime); + std::shared_ptr loadBalancer, // + double startTime, double endTime, amr::LoadBalancerDetails const& lb_info, + int loadBalancerPatchId); private: + amr::LoadBalancerDetails const lb_info_; + bool const is_tagging_refinement = false; + int loadBalancerPatchId_ = -1; + std::size_t rebalance_coarsest_auto_back_off = 0; + std::size_t time_step_idx = 0; + std::size_t rebalance_coarsest_auto_back_off_by = 1; + + std::shared_ptr timeRefIntegrator_; + + + // essentially CascadePartitioner::computeNonUniformWorkload which is private + auto computeNonUniformWorkLoadForLevel0() const + { + auto const& L0 = *timeRefIntegrator_->getPatchHierarchy()->getPatchLevel(0); + double load = 0.0; + for (SAMRAI::hier::PatchLevel::iterator ip(L0.begin()); ip != L0.end(); ++ip) + { + std::shared_ptr const& patch = *ip; + load += SAMRAI::mesh::BalanceUtilities::computeNonUniformWorkload( + patch, loadBalancerPatchId_, patch->getBox()); + } + return load; + } + + bool _should_rebalance_now(); + + bool cadence_rebalance_check() + { + return (time_step_idx == 0 and lb_info_.on_init) + or (lb_info_.every > 0 and time_step_idx > 0 + and time_step_idx % lb_info_.every == 0); + } + + bool tolerance_rebalance_check() + { + if (rebalance_coarsest_auto_back_off == 0) + { + PHARE_LOG_SCOPE(1, "Integrator::_should_rebalance_now::automatic"); + + auto workLoads = core::mpi::collect(computeNonUniformWorkLoadForLevel0()); + auto const max_value = *std::max_element(workLoads.begin(), workLoads.end()); + for (auto& workload : workLoads) + workload /= max_value; + auto const min_value = *std::min_element(workLoads.begin(), workLoads.end()); + assert(min_value <= 1); + + if ((1 - min_value) > lb_info_.tolerance) + { + rebalance_coarsest_auto_back_off_by = lb_info_.next_rebalance; + rebalance_coarsest_auto_back_off = rebalance_coarsest_auto_back_off_by; + return true; + } + + rebalance_coarsest_auto_back_off_by *= lb_info_.next_rebalance_backoff_multiplier; + rebalance_coarsest_auto_back_off = rebalance_coarsest_auto_back_off_by; + + if (rebalance_coarsest_auto_back_off > lb_info_.max_next_rebalance) + rebalance_coarsest_auto_back_off_by = rebalance_coarsest_auto_back_off + = lb_info_.max_next_rebalance; + } + + --rebalance_coarsest_auto_back_off; + return false; + }; + + std::function _rebalance_check; }; @@ -68,22 +153,23 @@ Integrator<_dimension>::Integrator( PHARE::initializer::PHAREDict const& dict, std::shared_ptr hierarchy, std::shared_ptr timeRefLevelStrategy, - std::shared_ptr tagAndInitStrategy, double startTime, - double endTime) + std::shared_ptr tagAndInitStrategy, + std::shared_ptr loadBalancer, double startTime, + double endTime, amr::LoadBalancerDetails const& lb_info, int loadBalancerPatchId) + : lb_info_{lb_info} + , is_tagging_refinement{_is_tagging_refinement(dict)} + , loadBalancerPatchId_{loadBalancerPatchId} + , rebalance_coarsest_auto_back_off{lb_info_.on_init ? 0 : lb_info_.next_rebalance} + , _rebalance_check{lb_info_.automatic ? std::bind(&Integrator::tolerance_rebalance_check, this) + : std::bind(&Integrator::cadence_rebalance_check, this)} { - // auto loadBalancer = std::make_shared( - // SAMRAI::tbox::Dimension{dimension}, "LoadBalancer"); - - auto loadBalancer = std::make_shared( - SAMRAI::tbox::Dimension{dimension}, "LoadBalancer"); - - loadBalancer->setSAMRAI_MPI(SAMRAI::tbox::SAMRAI_MPI::getSAMRAIWorld()); // TODO Is it really needed ? + loadBalancer->setSAMRAI_MPI( + SAMRAI::tbox::SAMRAI_MPI::getSAMRAIWorld()); // TODO Is it really needed ? auto refineDB = getUserRefinementBoxesDatabase(dict["simulation"]["AMR"]); auto standardTag = std::make_shared( "StandardTagAndInitialize", tagAndInitStrategy.get(), refineDB); - auto clustering = [&]() -> std::shared_ptr { if (!dict["simulation"]["AMR"].contains("clustering")) throw std::runtime_error(std::string{"clustering type not specificed"}); @@ -126,8 +212,14 @@ Integrator<_dimension>::Integrator( "TimeRefinementIntegrator", db, hierarchy, timeRefLevelStrategy, gridding); } +template +bool Integrator<_dimension>::_should_rebalance_now() +{ + if (is_tagging_refinement and lb_info_.active) + return _rebalance_check(); - + return false; +} template std::shared_ptr diff --git a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp index caf21df97..49098af06 100644 --- a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp +++ b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp @@ -42,13 +42,15 @@ class MaxwellianParticleInitializer : public ParticleInitializer const& bulkVelocity, std::array const& thermalVelocity, double const particleCharge, std::uint32_t const& nbrParticlesPerCell, std::optional seed = {}, - Basis const basis = Basis::Cartesian, - std::array const magneticField = {nullptr, nullptr, nullptr}) + Basis const basis = Basis::Cartesian, + std::array const& magneticField = {nullptr, nullptr, nullptr}, + double densityCutOff = 1e-5) : density_{density} , bulkVelocity_{bulkVelocity} , thermalVelocity_{thermalVelocity} , magneticField_{magneticField} , particleCharge_{particleCharge} + , densityCutOff_{densityCutOff} , nbrParticlePerCell_{nbrParticlesPerCell} , basis_{basis} , rngSeed_{seed} @@ -86,6 +88,7 @@ class MaxwellianParticleInitializer : public ParticleInitializer magneticField_; double particleCharge_; + double densityCutOff_ = 1e-5; std::uint32_t nbrParticlePerCell_; Basis basis_; std::optional rngSeed_; @@ -100,6 +103,7 @@ class MaxwellianInitFunctions FunctionArray& thermalVelocity, FunctionArray& magneticField, Basis const& basis, Coords const&... coords) : _n{density(coords...)} + { static_assert(sizeof...(coords) <= 3, "can only provide up to 3 coordinates"); for (std::uint32_t i = 0; i < 3; i++) @@ -123,7 +127,7 @@ class MaxwellianInitFunctions return {v[0]->data(), v[1]->data(), v[2]->data()}; } - std::shared_ptr> const _n; + std::shared_ptr> const _n, _ppc; std::array>, 3> _B, _V, _Vth; }; @@ -174,11 +178,14 @@ void MaxwellianParticleInitializer::loadParticles( auto randGen = getRNG(rngSeed_); ParticleDeltaDistribution deltaDistrib; - for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); flatCellIdx++) + for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); ++flatCellIdx) { + if (n[flatCellIdx] < densityCutOff_) + continue; + auto const cellWeight = n[flatCellIdx] / nbrParticlePerCell_; auto const AMRCellIndex = layout.localToAMR(point(flatCellIdx, ndCellIndices)); - + auto const iCell = AMRCellIndex.template toArray(); std::array particleVelocity; std::array, 3> basis; @@ -197,8 +204,7 @@ void MaxwellianParticleInitializer::loadParticles( if (basis_ == Basis::Magnetic) particleVelocity = basisTransform(basis, particleVelocity); - particles.emplace_back(Particle{cellWeight, particleCharge_, - AMRCellIndex.template toArray(), + particles.emplace_back(Particle{cellWeight, particleCharge_, iCell, deltas(deltaDistrib, randGen), particleVelocity}); } } diff --git a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp b/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp index 2f4cfa3f2..97cfa6b0c 100644 --- a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp +++ b/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp @@ -47,6 +47,8 @@ namespace core auto charge = dict["charge"].template to(); + auto densityCutOff = cppdict::get_value(dict, "density_cut_off", double{1e-16}); + auto nbrPartPerCell = dict["nbr_part_per_cell"].template to(); auto basisName = dict["basis"].template to(); @@ -59,22 +61,25 @@ namespace core if (dict.contains("init") && dict["init"].contains("seed")) seed = dict["init"]["seed"].template to>(); + std::array magneticField = {nullptr, nullptr, nullptr}; + if (basisName == "cartesian") { return std::make_unique< MaxwellianParticleInitializer>( - density, v, vth, charge, nbrPartPerCell, seed); + density, v, vth, charge, nbrPartPerCell, seed, Basis::Cartesian, + magneticField, densityCutOff); } else if (basisName == "magnetic") { - [[maybe_unused]] Basis basis = Basis::Magnetic; - [[maybe_unused]] auto& bx = dict["magnetic_x"].template to(); - [[maybe_unused]] auto& by = dict["magnetic_x"].template to(); - [[maybe_unused]] auto& bz = dict["magnetic_x"].template to(); + magneticField[0] = dict["magnetic_x"].template to(); + magneticField[1] = dict["magnetic_x"].template to(); + magneticField[2] = dict["magnetic_x"].template to(); return std::make_unique< MaxwellianParticleInitializer>( - density, v, vth, charge, nbrPartPerCell, seed); + density, v, vth, charge, nbrPartPerCell, seed, Basis::Magnetic, + magneticField, densityCutOff); } } // TODO throw? diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index c9bef9837..940fabe60 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -186,6 +186,17 @@ class NdArrayView : NdArrayViewer }; +template +auto make_array_view(DataType* data, std::array shape) +{ + return NdArrayView{data, shape}; +} + +template +auto make_array_view(DataType const* const data, std::array shape) +{ + return NdArrayView{data, shape}; +} template diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index fd3e1fe09..87472e75c 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -158,8 +158,8 @@ class box_iterator } public: - Point operator*() { return index_; } - + auto& operator*() const { return index_; } + auto operator->() const { return &index_; } void increment(std::size_t idim) { diff --git a/src/initializer/data_provider.hpp b/src/initializer/data_provider.hpp index 24ce7a0f5..c48f31a6e 100644 --- a/src/initializer/data_provider.hpp +++ b/src/initializer/data_provider.hpp @@ -48,9 +48,9 @@ namespace initializer using InitFunction = typename InitFunctionHelper::type; - using PHAREDict = cppdict::Dict, double, std::vector, std::size_t, - std::optional, std::string, InitFunction<1>, - InitFunction<2>, InitFunction<3>>; + using PHAREDict = cppdict::Dict, double, std::vector, + std::size_t, std::optional, std::string, + InitFunction<1>, InitFunction<2>, InitFunction<3>>; class PHAREDictHandler diff --git a/src/initializer/dictator.cpp b/src/initializer/dictator.cpp index 21c4d7b37..5d41fa7e8 100644 --- a/src/initializer/dictator.cpp +++ b/src/initializer/dictator.cpp @@ -45,6 +45,7 @@ PYBIND11_MODULE(dictator, m) m.def("add_size_t", add, "add_size_t"); m.def("add_optional_size_t", add>, "add_optional_size_t"); + m.def("add_bool", add, "add"); m.def("add_int", add, "add"); m.def("add_vector_int", add>, "add"); m.def("add_double", add, "add"); diff --git a/src/simulator/simulator.hpp b/src/simulator/simulator.hpp index 2ec834912..e13069185 100644 --- a/src/simulator/simulator.hpp +++ b/src/simulator/simulator.hpp @@ -13,10 +13,10 @@ #include "core/utilities/mpi_utils.hpp" #include "core/utilities/timestamps.hpp" #include "amr/tagging/tagger_factory.hpp" +#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" - namespace PHARE { class ISimulator @@ -70,7 +70,9 @@ class Simulator : public ISimulator } if (dMan) + { return dMan->dump(timestamp, timestep); + } return false; } @@ -261,27 +263,40 @@ void Simulator::hybrid_init(initializer::PHAREDict auto hybridTagger_ = amr::TaggerFactory::make("HybridModel", "default"); multiphysInteg_->registerTagger(0, maxLevelNumber_ - 1, std::move(hybridTagger_)); - - + amr::LoadBalancerDetails lb_info + = amr::LoadBalancerDetails::FROM(dict["simulation"]["AMR"]["loadbalancing"]); auto lbm_ = std::make_unique>(dict); + auto lbe_ = std::make_shared>(lb_info.mode, + lbm_->getId()); + + auto loadBalancer_db = std::make_shared("LoadBalancerDB"); + loadBalancer_db->putDouble("flexible_load_tolerance", lb_info.tolerance); + auto loadBalancer = std::make_shared( + SAMRAI::tbox::Dimension{dimension}, "LoadBalancer", loadBalancer_db); + + if (dict["simulation"]["AMR"]["refinement"].contains("tagging")) + { // Load balancers break with refinement boxes - only tagging supported + /* + P=0000000:Program abort called in file ``/.../SAMRAI/xfer/RefineSchedule.cpp'' at line 369 + P=0000000:ERROR MESSAGE: + P=0000000:RefineSchedule:RefineSchedule error: We are not currently + P=0000000:supporting RefineSchedules with the source level finer + P=0000000:than the destination level + */ + lbm_->addLoadBalancerEstimator(0, maxLevelNumber_ - 1, std::move(lbe_)); + lbm_->setLoadBalancer(loadBalancer); + } - auto lbe_ = std::make_unique>( - dict["simulation"]["AMR"]["loadbalancing"].template to(), lbm_->getId()); - - lbm_->addLoadBalancerEstimator(0, maxLevelNumber_ - 1, std::move(lbe_)); - lbm_->addLoadBalancer(std::make_unique( - SAMRAI::tbox::Dimension{dim}, "cascade")); + auto lbm_id = lbm_->getId(); // moved on next line multiphysInteg_->setLoadBalancerManager(std::move(lbm_)); - - - if (dict["simulation"].contains("restarts")) startTime_ = restarts_init(dict["simulation"]["restarts"]); - integrator_ = std::make_unique(dict, hierarchy_, multiphysInteg_, multiphysInteg_, - startTime_, finalTime_); + integrator_ + = std::make_unique(dict, hierarchy_, multiphysInteg_, multiphysInteg_, + loadBalancer, startTime_, finalTime_, lb_info, lbm_id); timeStamper = core::TimeStamperFactory::create(dict["simulation"]); @@ -374,6 +389,7 @@ void Simulator<_dimension, _interp_order, _nbRefinedPart>::initialize() template double Simulator<_dimension, _interp_order, _nbRefinedPart>::advance(double dt) { + PHARE_LOG_SCOPE(1, "Simulator::advance"); double dt_new = 0; if (!integrator_) @@ -382,6 +398,7 @@ double Simulator<_dimension, _interp_order, _nbRefinedPart>::advance(double dt) try { PHARE_LOG_SCOPE(1, "Simulator::advance"); + dt_new = integrator_->advance(dt); currentTime_ = startTime_ + ((*timeStamper) += dt); } diff --git a/tests/core/data/ions/test_ions.cpp b/tests/core/data/ions/test_ions.cpp index 77192a467..19d7e62c5 100644 --- a/tests/core/data/ions/test_ions.cpp +++ b/tests/core/data/ions/test_ions.cpp @@ -103,7 +103,6 @@ class theIons : public ::testing::Test dict["ions"]["pop1"]["particle_initializer"]["thermal_velocity_z"] = static_cast(vthz); - dict["ions"]["pop1"]["particle_initializer"]["nbrPartPerCell"] = int{100}; dict["ions"]["pop1"]["particle_initializer"]["charge"] = -1.; dict["ions"]["pop1"]["particle_initializer"]["basis"] = std::string{"Cartesian"}; diff --git a/tests/diagnostic/__init__.py b/tests/diagnostic/__init__.py index 24fddb61d..a4c4b0cb1 100644 --- a/tests/diagnostic/__init__.py +++ b/tests/diagnostic/__init__.py @@ -20,6 +20,12 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None): if timestamps is None: timestamps = all_timestamps(sim) + ph.InfoDiagnostics( + quantity="particle_count", + write_timestamps=timestamps, + compute_timestamps=timestamps, + ) + ph.MetaDiagnostics( quantity="tags", write_timestamps=timestamps, diff --git a/tests/functional/harris/harris_2d_lb.py b/tests/functional/harris/harris_2d_lb.py new file mode 100644 index 000000000..e701e8231 --- /dev/null +++ b/tests/functional/harris/harris_2d_lb.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 + +import os +import numpy as np +import matplotlib as mpl +from pathlib import Path + +import pyphare.pharein as ph +from pyphare.cpp import cpp_lib +from pyphare.pharesee.run import Run +from pyphare.simulator.simulator import Simulator, startMPI + +from tests.simulator import SimulatorTest +from tools.python3 import plotting as m_plotting + +mpl.use("Agg") + +LOAD_BALANCE = os.getenv("LOAD_BALANCE", "True").lower() in ("true", "1", "t") + +cpp = cpp_lib() +startMPI() + +cells = (800, 800) +time_step = 0.005 +final_time = 50 +timestamps = np.arange(0, final_time + time_step, final_time / 5) + +if cpp.mpi_rank() == 0: + print(LOAD_BALANCE, "diag timestamps:", timestamps) + +diag_dir = "phare_outputs/harris_lb" +if not LOAD_BALANCE: + diag_dir = "phare_outputs/harris" + +plot_dir = Path(f"{diag_dir}_plots") +plot_dir.mkdir(parents=True, exist_ok=True) + + +def config(): + L = 0.5 + + sim = ph.Simulation( + time_step=time_step, + final_time=final_time, + cells=cells, + dl=(0.40, 0.40), + refinement="tagging", + max_nbr_levels=2, + nesting_buffer=1, + clustering="tile", + tag_buffer="1", + hyper_resistivity=0.002, + resistivity=0.001, + diag_options={ + "format": "phareh5", + "options": {"dir": diag_dir, "mode": "overwrite"}, + }, + restart_options={ + "dir": "checkpoints", + "mode": "overwrite", + "timestamps": timestamps, + # "restart_time": 0.0, + }, + ) + + def density(x, y): + Ly = sim.simulation_domain()[1] + return ( + 0.4 + + 1.0 / np.cosh((y - Ly * 0.3) / L) ** 2 + + 1.0 / np.cosh((y - Ly * 0.7) / L) ** 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] + sigma = 1.0 + dB = 0.1 + + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + + dBy1 = 2 * dB * x0 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2) + dBy2 = -2 * dB * x0 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2) + + return dBy1 + dBy2 + + def bx(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + sigma = 1.0 + dB = 0.1 + + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + + dBx1 = -2 * dB * y1 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2) + dBx2 = 2 * dB * y2 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2) + + v1 = -1 + v2 = 1.0 + return v1 + (v2 - v1) * (S(y, Ly * 0.3, L) - S(y, Ly * 0.7, L)) + dBx1 + dBx2 + + 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 = 0.7 + temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5) + assert np.all(temp > 0) + return temp + + def vxyz(x, y): + return 0.0 + + def vthxyz(x, y): + return np.sqrt(T(x, y)) + + vvv = {**{f"vbulk{c}": vxyz for c in "xyz"}, **{f"vth{c}": vthxyz for c in "xyz"}} + + ph.MaxwellianFluidModel( + bx=bx, by=by, bz=bz, protons={"charge": 1, "density": density, **vvv} + ) + ph.ElectronModel(closure="isothermal", Te=0.0) + + for quantity in ["E", "B"]: + ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) + for quantity in ["density", "bulkVelocity"]: + ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) + + ph.FluidDiagnostics( + quantity="density", write_timestamps=timestamps, population_name="protons" + ) + ph.InfoDiagnostics(quantity="particle_count") + + if LOAD_BALANCE: + ph.LoadBalancer(active=True, auto=True, mode="nppc", tol=0.05) + + return sim + + +def plot_file_for_qty(qty, time): + return f"{plot_dir}/harris_{qty}_t{time}.png" + + +def plot(diag_dir): + run = Run(diag_dir) + for time in timestamps: + run.GetDivB(time).plot( + filename=plot_file_for_qty("divb", time), + plot_patches=True, + vmin=1e-11, + vmax=2e-10, + ) + run.GetRanks(time).plot( + filename=plot_file_for_qty("Ranks", time), + plot_patches=True, + ) + run.GetN(time, pop_name="protons").plot( + filename=plot_file_for_qty("N", time), + plot_patches=True, + ) + for c in ["x", "y", "z"]: + run.GetB(time).plot( + filename=plot_file_for_qty(f"b{c}", time), + qty=f"B{c}", + plot_patches=True, + ) + run.GetJ(time).plot( + filename=plot_file_for_qty("jz", time), + qty="Jz", + plot_patches=True, + vmin=-2, + vmax=2, + ) + + +class HarrisTest(SimulatorTest): + def __init__(self, *args, **kwargs): + super(HarrisTest, self).__init__(*args, **kwargs) + self.simulator = None + + def tearDown(self): + super(HarrisTest, self).tearDown() + if self.simulator is not None: + self.simulator.reset() + self.simulator = None + ph.global_vars.sim = None + + def test_run(self): + self.register_diag_dir_for_cleanup(diag_dir) + Simulator(config()).run().reset() + if cpp.mpi_rank() == 0: + plot(diag_dir) + m_plotting.plot_run_timer_data(diag_dir, cpp.mpi_rank()) + cpp.mpi_barrier() + return self + + +if __name__ == "__main__": + HarrisTest().test_run().tearDown() diff --git a/tests/simulator/CMakeLists.txt b/tests/simulator/CMakeLists.txt index 09aecb003..86274d16f 100644 --- a/tests/simulator/CMakeLists.txt +++ b/tests/simulator/CMakeLists.txt @@ -22,6 +22,9 @@ if(HighFive) if(testMPI) phare_mpi_python3_exec(9 3 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) phare_mpi_python3_exec(9 4 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) + + # doesn't make sense in serial + phare_mpi_python3_exec(9 3 load_balancing test_load_balancing.py ${CMAKE_CURRENT_BINARY_DIR}) endif(testMPI) phare_python3_exec(11, test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index 56ca9e9d9..11ea9f8e1 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -1,3 +1,4 @@ +import os import unittest from datetime import datetime import pyphare.pharein as ph, numpy as np @@ -191,6 +192,55 @@ def _diff(slice0): return boxes +# +# + + +def caliper_func_times_json(data_dir, mpi_rank=0): + return f"{os.path.join(data_dir, f'func_times.{mpi_rank}.json')}" + + +def caliper_recorder_cali(data_dir, mpi_rank=0): + return f"{os.path.join(data_dir, f'recorder.{mpi_rank}.cali')}" + + +CALIPER_MODES = [ + # "callpath:event:recorder:trace", + "report,event,trace,timestamp,recorder", # light + "alloc,aggregate,cpuinfo,memusage,debug,env,event,loop_monitor,region_monitor,textlog,io,pthread,sysalloc,recorder,report,timestamp,statistics,spot,trace,validator,mpi,mpireport,mpiflush", # heavy +] + + +def activate_caliper(data_dir, mode_idx=0): + from pyphare.cpp import cpp_lib + + rank = cpp_lib().mpi_rank() + env = os.environ + + # env["CALI_SERVICES_ENABLE"] = "event,trace,timer,report" + # env["CALI_REPORT_CONFIG"] = "format json" + # env["CALI_REPORT_FILENAME"] = "trace.json" + + # env[ + # "CALI_CONFIG" + # ] = "hatchet-region-profile,topdown-counters.all,output.format=json" + + # # env["CALI_CONFIG_PROFILE"] = "callstack-trace" + # env["CALI_SERVICES_ENABLE"] = CALIPER_MODES[mode_idx] + + # env["CALI_CONFIG"] = "hatchet-region-profile" + + # # env["CALI_CALLPATH_USE_NAME"] = "true" + + # env["CALI_REPORT_FILENAME"] = caliper_func_times_json(data_dir, rank) + # env[ + # "CALI_REPORT_CONFIG" + # ] = "SELECT function,time.duration ORDER BY time.duration FORMAT json" + # env["CALI_RECORDER_FILENAME"] = caliper_recorder_cali(data_dir, rank) + + # print("os.environ", os.environ) + + class SimulatorTest(unittest.TestCase): test_kwargs = ["rethrow"] diff --git a/tests/simulator/test_load_balancing.py b/tests/simulator/test_load_balancing.py new file mode 100644 index 000000000..1dbb89221 --- /dev/null +++ b/tests/simulator/test_load_balancing.py @@ -0,0 +1,284 @@ +#!/usr/bin/env python3 + +# basically a harris run + +import unittest +import pyphare.pharein as ph +from pyphare.pharesee.hierarchy import hierarchy_compare +from pyphare.pharesee.particles import single_patch_per_level_per_pop_from + + +from pyphare.simulator.simulator import Simulator, startMPI +from tests.simulator import SimulatorTest + +import numpy as np +import matplotlib as mpl + +mpl.use("Agg") + +from pyphare.cpp import cpp_lib + +cpp = cpp_lib() +startMPI() + +ndim = 2 +interp = 1 +mpi_size = cpp.mpi_size() +time_step_nbr = 2 +time_step = 0.001 +cells = (100, 100) +dl = (0.2, 0.2) +diag_outputs = "phare_outputs/load_balancing_2d" +timestamps = [x * time_step for x in range(time_step_nbr + 1)] + + +def config(diag_dir, loadbalancing={}): + sim = ph.Simulation( + time_step_nbr=time_step_nbr, + time_step=time_step, + cells=cells, + dl=dl, + refinement="tagging", + max_nbr_levels=2, + hyper_resistivity=0.001, + resistivity=0.001, + diag_options={ + "format": "phareh5", + "options": {"dir": diag_dir, "mode": "overwrite"}, + }, + ) + + 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 densityHigh(x, y): + assert cells == (100, 100) # next lines are dependent + rho = y.copy() + rho[:] = 1e-10 + rho[np.where(np.isclose(y, 10, atol=0.1))] = 1 + return rho + + def densityLow(x, y): + assert cells == (100, 100) # next lines are dependent + rho = y.copy() + rho[:] = 0.1 + rho[np.where(np.isclose(y, 10, atol=0.1))] = 0 + return rho + + 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 vxyz(x, y): + return 0.0 + + def vthxyz(x, y): + return np.sqrt(T(x, y)) + + vHigh = { + "vbulkx": vxyz, + "vbulky": vxyz, + "vbulkz": vxyz, + "vthx": vthxyz, + "vthy": vthxyz, + "vthz": vthxyz, + "nbr_part_per_cell": 590, + "density_cut_off": 1e-5, + } + + vLow = { + "vbulkx": vxyz, + "vbulky": vxyz, + "vbulkz": vxyz, + "vthx": vthxyz, + "vthy": vthxyz, + "vthz": vthxyz, + "nbr_part_per_cell": 90, + "density_cut_off": 1e-5, + } + + ph.MaxwellianFluidModel( + bx=bx, + by=by, + bz=bz, + brotons={"charge": 1, "density": densityHigh, **vHigh, "init": {"seed": 12334}}, + protons={"charge": 1, "density": densityLow, **vLow, "init": {"seed": 43210}}, + ) + ph.ElectronModel(closure="isothermal", Te=0.0) + + for pop in ["brotons", "protons"]: + ph.ParticleDiagnostics( + quantity="domain", write_timestamps=timestamps, population_name=pop + ) + + if loadbalancing: + ph.LoadBalancer(**loadbalancing) + + return sim + + +def get_time(path, time=0, pop="protons", datahier=None): + time = "{:.10f}".format(time) + from pyphare.pharesee.hierarchy import hierarchy_from + + return hierarchy_from( + h5_filename=path + f"/ions_pop_{pop}_domain.h5", time=time, hier=datahier + ) + + +def get_particles(diag_dir, time=0): + hier = get_time(diag_dir, time) + hier = get_time(diag_dir, time, "brotons", hier) + return hier + + +def time_info(diag_dir, time=0): + hier = get_particles(diag_dir, time) + + per_rank = {f"p{rank}": 0 for rank in range(mpi_size)} + + def _parse_rank(patch_id): + return patch_id.split("#")[0] + + for ilvl, lvl in hier.levels().items(): + for patch in lvl: + for pd_key, pd in patch.patch_datas.items(): + per_rank[_parse_rank(patch.id)] += pd.size() + + return per_rank + + +class LoadBalancingTest(SimulatorTest): + def tearDown(self): + ph.global_vars.sim = None + + def run_sim(self, diags_dir, dic={}): + ph.global_vars.sim = None + self.register_diag_dir_for_cleanup(diags_dir) + Simulator(config(diags_dir, dic)).run() + return diags_dir + + def test_has_balanced_auto(self): + if mpi_size == 1: # doesn't make sense + return + + diag_dir = self.run_sim( + self.unique_diag_dir_for_test_case(diag_outputs, ndim, interp), + dict(active=True, auto=True, mode="nppc", tol=0.05), + ) + + if cpp.mpi_rank() > 0: + return + + t0_sdev = np.std(list(time_info(diag_dir).values())) + tend_sdev = np.std(list(time_info(diag_dir, timestamps[-1]).values())) + self.assertLess(tend_sdev, t0_sdev * 0.1) # empirical + + def test_has_balanced_manual(self): + if mpi_size == 1: # doesn't make sense + return + + diag_dir = self.run_sim( + self.unique_diag_dir_for_test_case(diag_outputs, ndim, interp), + dict(active=True, on_init=True, every=1, mode="nppc", tol=0.05), + ) + + if cpp.mpi_rank() > 0: + return + + t0_sdev = np.std(list(time_info(diag_dir).values())) + tend_sdev = np.std(list(time_info(diag_dir, timestamps[-1]).values())) + self.assertLess(tend_sdev, t0_sdev * 0.1) # empirical + + def test_has_not_balanced_as_defaults(self): + if mpi_size == 1: # doesn't make sense + return + + diag_dir = self.run_sim( + self.unique_diag_dir_for_test_case(diag_outputs, ndim, interp) + ) + + if cpp.mpi_rank() > 0: + return + + t0_sdev = np.std(list(time_info(diag_dir).values())) + tend_sdev = np.std(list(time_info(diag_dir, timestamps[-1]).values())) + self.assertGreater(tend_sdev, t0_sdev * 0.1) # empirical + + def test_compare_is_and_is_not_balanced(self): + if mpi_size == 1: # doesn't make sense + return + + check_time = 0.001 + + not_hier = get_particles( + self.run_sim( + self.unique_diag_dir_for_test_case(diag_outputs + "/not", ndim, interp) + ), + check_time, + ) + is_hier = get_particles( + self.run_sim( + self.unique_diag_dir_for_test_case(diag_outputs, ndim, interp), + dict(active=True, auto=True, mode="nppc", tol=0.05), + ), + check_time, + ) + + if cpp.mpi_rank() == 0: + not_hier = single_patch_per_level_per_pop_from(not_hier) + is_hier = single_patch_per_level_per_pop_from(is_hier) + self.assertTrue(hierarchy_compare(not_hier, is_hier)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/simulator/test_tagging.py b/tests/simulator/test_tagging.py.off similarity index 99% rename from tests/simulator/test_tagging.py rename to tests/simulator/test_tagging.py.off index 11c3af280..0785dcf77 100644 --- a/tests/simulator/test_tagging.py +++ b/tests/simulator/test_tagging.py.off @@ -1,5 +1,7 @@ #!/usr/bin/env python3 +# renamed to not be run on test discovery by directory + from pyphare.cpp import cpp_lib cpp = cpp_lib() diff --git a/tools/python3/__init__.py b/tools/python3/__init__.py index 16db65902..fd321517a 100644 --- a/tools/python3/__init__.py +++ b/tools/python3/__init__.py @@ -11,11 +11,14 @@ def run(cmd, shell=True, capture_output=True, check=False, print_cmd=True, **kwa if print_cmd: print(f"running: {cmd}") try: - return subprocess.run(cmd, shell=shell, capture_output=capture_output, check=check, **kwargs) - except subprocess.CalledProcessError as e: # only triggers on failure if check=True + return subprocess.run( + cmd, shell=shell, capture_output=capture_output, check=check, **kwargs + ) + except subprocess.CalledProcessError as e: # only triggers on failure if check=True print(f"run failed with error: {e}\n\t{e.stdout}\n\t{e.stderr} ") raise RuntimeError(decode_bytes(e.stderr)) + def run_mp(cmds, N_CORES=None, **kwargs): """ spawns N_CORES threads (default=len(cmds)) running commands and waiting for results @@ -65,13 +68,46 @@ def scan_dir(path, files_only=False, dirs_only=False, drop=[]): if all([check(entry) for check in checks]) ] + import contextlib + + @contextlib.contextmanager def pushd(new_cwd): import os + cwd = os.getcwd() os.chdir(new_cwd) try: yield finally: os.chdir(cwd) + + +@contextlib.contextmanager +def extend_env_path(paths): + import os + + if isinstance(paths, str): + paths = [paths] + old_path = os.environ["PATH"] + for path in paths: + os.environ["PATH"] += os.pathsep + path + try: + yield + finally: + os.environ["PATH"] = old_path + + +@contextlib.contextmanager +def extend_sys_path(paths): + import sys + + if isinstance(paths, str): + paths = [paths] + old_path = sys.path[:] + sys.path.extend(paths) + try: + yield + finally: + sys.path = old_path diff --git a/tools/python3/caliper.py b/tools/python3/caliper.py new file mode 100644 index 000000000..86e423d14 --- /dev/null +++ b/tools/python3/caliper.py @@ -0,0 +1,149 @@ +# +# +# + +from pathlib import Path +from tools.python3 import extend_env_path, run + +SCRIPT_DIR = Path(__file__).resolve().parent + +# /home/deegan/git/phare/master/build/subprojects/caliper/src/tools/cali-query +PROJECT_ROOT_DIR = SCRIPT_DIR.parent.parent +DEFAULT_CALI_QUERY_PATH = ( + f"{PROJECT_ROOT_DIR}/build/subprojects/caliper/src/tools/cali-query" +) + + +# https://hatchet.readthedocs.io/en/latest/basic_tutorial.html#introduction +def hatchet_on_caliper(cali_file="phare_outputs/recorder.0.cali"): + import hatchet as ht + import subprocess + + cali_file = "region_profile.cali" + + # Setup desired cali query. + # cali_query = "cali-query" + # grouping_attribute = "function" + # default_metric = "sum(sum#time.duration),inclusive_sum(sum#time.duration)" + # query = "select function,%s group by %s format json-split" % ( + # default_metric, + # grouping_attribute, + # ) + + # # Use ``cali-query`` here to produce the json-split stream. + # with extend_env_path(DEFAULT_CALI_QUERY_PATH): + # cali_json = subprocess.Popen( + # [cali_query, "-q", query, cali_file], stdout=subprocess.PIPE + # ) + + # # Use hatchet's ``from_caliper`` API with the resulting json-split. + # # The result is stored into Hatchet's GraphFrame. + # gf = ht.GraphFrame.from_caliper(cali_json.stdout) + + # # Printout the DataFrame component of the GraphFrame. + # print(gf.dataframe) + + # # Printout the graph component of the GraphFrame. + # # Use "time (inc)" as the metric column to be displayed + # print(gf.tree(metric_column="time (inc)")) + + # # Setup desired cali query. + # grouping_attribute = "function" + # default_metric = "sum(sum#time.duration),inclusive_sum(sum#time.duration)" + # query = "select function,%s group by %s format json-split" % ( + # default_metric, + # grouping_attribute, + # ) + + # # Use hatchet's ``from_caliper`` API with the path to the cali file and the + # # query. This API will internally run ``cali-query`` on this file to + # # produce a json-split stream. The result is stored into Hatchet's + # # GraphFrame. + # with extend_env_path(DEFAULT_CALI_QUERY_PATH): + # gf = ht.GraphFrame.from_caliper(cali_file, query) + + # # Printout the DataFrame component of the GraphFrame. + # print(gf.dataframe) + + # # Printout the graph component of the GraphFrame. + # # Use "time (inc)" as the metric column to be displayed + # print(gf.tree(metric_column="time (inc)")) + + grouping_attribute = "function" + default_metric = "sum(sum#time.duration),inclusive_sum(sum#time.duration)" + query = "select function,%s group by %s format json-split" % ( + default_metric, + grouping_attribute, + ) + + query = "SELECT function,time.duration ORDER BY time.duration FORMAT json-split" + with extend_env_path(DEFAULT_CALI_QUERY_PATH): + gf = ht.GraphFrame.from_caliper(cali_file, query) + print(gf.dataframe) + + +def hatchet_region_profile_inclusive(time_per_fn_per_rank, ranks): + fn_ts = {r: {} for r in range(ranks)} + + for rank in range(ranks): + for k, v in time_per_fn_per_rank[rank].items(): + # print(k, v) + if "/" not in k: # is root scope + fn_ts[rank][k] = v + else: + bits = k.split("/") + fn_ts[rank][bits[0]] += v + + return fn_ts + + +def hatchet_region_profile(json_file="trace.json", ranks=2, inclusive=False): + import json + + time_per_rank = {r: 0 for r in range(ranks)} + time_per_fn_per_rank = {r: {} for r in range(ranks)} + time_per_fn_per_ts_per_rank = {r: [] for r in range(ranks)} + + ts_idx = 0 + post_init = False + with open(json_file, "r") as f: + for d in json.load(f): + if "event.end#region" not in d: + continue + rank = d["mpi.rank"] + time = d["time"] + time_per_rank[rank] += time + if "path" not in d: + continue + fn = d["path"] + + if post_init and fn == "Simulator::advance": + print("post_init") + ts_idx += 1 + elif not post_init and fn == "Simulator::advance": + post_init = True + + if fn not in time_per_fn_per_rank[rank]: + time_per_fn_per_rank[rank][fn] = 0 + time_per_fn_per_rank[rank][fn] += time + + if post_init: + if ts_idx == len(time_per_fn_per_ts_per_rank[rank]): + time_per_fn_per_ts_per_rank[rank].append({}) + time_per_fn_per_ts_per_rank[rank][ts_idx][fn] = time + + if inclusive: + time_per_fn_per_rank = hatchet_region_profile_inclusive( + time_per_fn_per_rank, ranks + ) + + print(len(time_per_fn_per_ts_per_rank[0])) + # for i, d in enumerate(time_per_fn_per_ts_per_rank[0]): + # for k, v in d.items(): + # print(i, k, v) + + print("time_per_rank", time_per_rank) + + +if __name__ == "__main__": + hatchet_on_caliper()