diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 6b7356b28..3cbc770a9 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -1,9 +1,7 @@ import os import sys import subprocess -import numpy as np -from pyphare.core.phare_utilities import is_scalar from .uniform_model import UniformModel from .maxwellian_fluid_model import MaxwellianFluidModel from .electron_model import ElectronModel @@ -14,11 +12,7 @@ MetaDiagnostics, InfoDiagnostics, ) -from .simulation import ( - Simulation, - serialize as serialize_sim, - deserialize as deserialize_sim, -) +from .simulation import Simulation from .load_balancer import LoadBalancer __all__ = [ @@ -31,6 +25,7 @@ "MetaDiagnostics", "InfoDiagnostics", "Simulation", + "LoadBalancer", ] # This exists to allow a condition variable for when we are running PHARE from C++ via phare-exe @@ -64,58 +59,6 @@ def NO_GUI(): mpl.use("Agg") -def getSimulation(): - from .global_vars import sim - - return sim - - -def _patch_data_ids(restart_file_dir): - """ - for restarts we save samrai patch data ids to the restart files, which we access from here - to tell samrai which patch datas to load from the restart file on restart - """ - from pyphare.cpp import cpp_etc_lib - - return cpp_etc_lib().patch_data_ids(restart_file_dir) - - -def _serialized_simulation_string(restart_file_dir): - from pyphare.cpp import cpp_etc_lib - - return cpp_etc_lib().serialized_simulation_string(restart_file_dir) - - -# converts scalars to array of expected size -# converts lists to arrays -class py_fn_wrapper: - def __init__(self, fn): - self.fn = fn - - def __call__(self, *xyz): - args = [np.asarray(arg) for arg in xyz] - ret = self.fn(*args) - if isinstance(ret, list): - ret = np.asarray(ret) - if is_scalar(ret): - ret = np.full(len(args[-1]), ret) - return ret - - -# Wrap calls to user init functions to turn C++ vectors to ndarrays, -# and returned ndarrays to C++ span -class fn_wrapper(py_fn_wrapper): - def __init__(self, fn): - super().__init__(fn) - - def __call__(self, *xyz): - from pyphare.cpp import cpp_etc_lib - - # convert numpy array to C++ SubSpan - # couples vector init functions to C++ - return cpp_etc_lib().makePyArrayWrapper(super().__call__(*xyz)) - - def clearDict(): """ dict may contain dangling references from a previous simulation unless cleared @@ -126,289 +69,12 @@ def clearDict(): def populateDict(): - from .global_vars import sim as simulation - import pybindlibs.dictator as pp - - # pybind complains if receiving wrong type - 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)) - - def add_size_t(path, val): - casted = int(val) - if casted < 0: - raise RuntimeError("pyphare.__init__::add_size_t received negative value") - pp.add_size_t(path, casted) - - def add_vector_int(path, val): - pp.add_vector_int(path, list(val)) - - add_string = pp.add_string - addInitFunction = getattr(pp, "addInitFunction{:d}".format(simulation.ndim) + "D") - - add_string("simulation/name", "simulation_test") - add_int("simulation/dimension", simulation.ndim) - - if simulation.smallest_patch_size is not None: - add_vector_int( - "simulation/AMR/smallest_patch_size", simulation.smallest_patch_size - ) - if simulation.largest_patch_size is not None: - add_vector_int( - "simulation/AMR/largest_patch_size", simulation.largest_patch_size - ) - - add_string("simulation/grid/layout_type", simulation.layout) - add_int("simulation/grid/nbr_cells/x", simulation.cells[0]) - add_double("simulation/grid/meshsize/x", simulation.dl[0]) - add_double("simulation/grid/origin/x", simulation.origin[0]) - add_string("simulation/grid/boundary_type/x", simulation.boundary_types[0]) - - if simulation.ndim > 1: - add_int("simulation/grid/nbr_cells/y", simulation.cells[1]) - add_double("simulation/grid/meshsize/y", simulation.dl[1]) - add_double("simulation/grid/origin/y", simulation.origin[1]) - add_string("simulation/grid/boundary_type/y", simulation.boundary_types[1]) - - if simulation.ndim > 2: - add_int("simulation/grid/nbr_cells/z", simulation.cells[2]) - add_double("simulation/grid/meshsize/z", simulation.dl[2]) - add_double("simulation/grid/origin/z", simulation.origin[2]) - add_string("simulation/grid/boundary_type/z", simulation.boundary_types[2]) - - add_int("simulation/interp_order", simulation.interp_order) - add_int("simulation/refined_particle_nbr", simulation.refined_particle_nbr) - add_double("simulation/time_step", simulation.time_step) - add_int("simulation/time_step_nbr", simulation.time_step_nbr) - - add_string("simulation/AMR/clustering", simulation.clustering) - add_int("simulation/AMR/max_nbr_levels", simulation.max_nbr_levels) - add_vector_int("simulation/AMR/nesting_buffer", simulation.nesting_buffer) - - add_int("simulation/AMR/tag_buffer", simulation.tag_buffer) - - refinement_boxes = simulation.refinement_boxes - - def as_paths(rb): - add_int("simulation/AMR/refinement/boxes/nbr_levels/", len(rb.keys())) - for level, boxes in rb.items(): - level_path = "simulation/AMR/refinement/boxes/" + level + "/" - add_int(level_path + "nbr_boxes/", int(len(boxes))) - for box_i, box in enumerate(boxes): - box_id = "B" + str(box_i) - lower = box.lower - upper = box.upper - box_lower_path_x = box_id + "/lower/x/" - box_upper_path_x = box_id + "/upper/x/" - add_int(level_path + box_lower_path_x, lower[0]) - add_int(level_path + box_upper_path_x, upper[0]) - if len(lower) >= 2: - box_lower_path_y = box_id + "/lower/y/" - box_upper_path_y = box_id + "/upper/y/" - add_int(level_path + box_lower_path_y, lower[1]) - add_int(level_path + box_upper_path_y, upper[1]) - if len(lower) == 3: - box_lower_path_z = box_id + "/lower/z/" - box_upper_path_z = box_id + "/upper/z/" - add_int(level_path + box_lower_path_z, lower[2]) - add_int(level_path + box_upper_path_z, upper[2]) - - if refinement_boxes is not None and simulation.refinement == "boxes": - as_paths(refinement_boxes) - elif simulation.refinement == "tagging": - add_string("simulation/AMR/refinement/tagging/method", "auto") - # the two following params are hard-coded for now - # they will become configurable when we have multi-models or several methods - # per model - add_string("simulation/AMR/refinement/tagging/model", "HybridModel") - add_string("simulation/AMR/refinement/tagging/method", "default") - add_double( - "simulation/AMR/refinement/tagging/threshold", simulation.tagging_threshold - ) - else: - add_string( - "simulation/AMR/refinement/tagging/method", "none" - ) # integrator.h might want some looking at - - add_string("simulation/algo/ion_updater/pusher/name", simulation.particle_pusher) - - add_double("simulation/algo/ohm/resistivity", simulation.resistivity) - add_double("simulation/algo/ohm/hyper_resistivity", simulation.hyper_resistivity) - add_string("simulation/algo/ohm/hyper_mode", simulation.hyper_mode) - - # load balancer block start - lb = simulation.load_balancer or LoadBalancer(active=False, _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 - - if init_model.nbr_populations() < 0: - raise RuntimeError("Number of populations cannot be negative") - add_size_t("simulation/ions/nbrPopulations", init_model.nbr_populations()) - - partinit = "particle_initializer" - for pop_index, pop in enumerate(init_model.populations): - pop_path = "simulation/ions/pop" - partinit_path = pop_path + "{:d}/".format(pop_index) + partinit + "/" - d = modelDict[pop] - add_string(pop_path + "{:d}/name".format(pop_index), pop) - add_double(pop_path + "{:d}/mass".format(pop_index), d["mass"]) - add_string(partinit_path + "name", "maxwellian") - - addInitFunction(partinit_path + "density", fn_wrapper(d["density"])) - addInitFunction(partinit_path + "bulk_velocity_x", fn_wrapper(d["vx"])) - addInitFunction(partinit_path + "bulk_velocity_y", fn_wrapper(d["vy"])) - addInitFunction(partinit_path + "bulk_velocity_z", fn_wrapper(d["vz"])) - 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_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") - - add_string("simulation/electromag/magnetic/name", "B") - maginit_path = "simulation/electromag/magnetic/initializer/" - addInitFunction(maginit_path + "x_component", fn_wrapper(modelDict["bx"])) - addInitFunction(maginit_path + "y_component", fn_wrapper(modelDict["by"])) - addInitFunction(maginit_path + "z_component", fn_wrapper(modelDict["bz"])) - - serialized_sim = serialize_sim(simulation) - - #### adding diagnostics - - diag_path = "simulation/diagnostics/" - for diag in list(simulation.diagnostics.values()): - diag.attributes["serialized_simulation"] = serialized_sim - - type_path = diag_path + diag.type + "/" - name_path = type_path + diag.name - add_string(name_path + "/" + "type", diag.type) - add_string(name_path + "/" + "quantity", diag.quantity) - add_size_t(name_path + "/" + "flush_every", diag.flush_every) - pp.add_array_as_vector( - name_path + "/" + "write_timestamps", diag.write_timestamps - ) - pp.add_array_as_vector( - name_path + "/" + "compute_timestamps", diag.compute_timestamps - ) - - add_size_t(name_path + "/" + "n_attributes", len(diag.attributes)) - for attr_idx, attr_key in enumerate(diag.attributes): - add_string(name_path + "/" + f"attribute_{attr_idx}_key", attr_key) - add_string( - name_path + "/" + f"attribute_{attr_idx}_value", - diag.attributes[attr_key], - ) - - if len(simulation.diagnostics) > 0: - if simulation.diag_options is not None and "options" in simulation.diag_options: - add_string( - diag_path + "filePath", simulation.diag_options["options"]["dir"] - ) - if "mode" in simulation.diag_options["options"]: - add_string( - diag_path + "mode", simulation.diag_options["options"]["mode"] - ) - if "fine_dump_lvl_max" in simulation.diag_options["options"]: - add_int( - diag_path + "fine_dump_lvl_max", - simulation.diag_options["options"]["fine_dump_lvl_max"], - ) - else: - add_string(diag_path + "filePath", "phare_output") - #### diagnostics added - - #### adding restarts - if simulation.restart_options is not None: - restart_options = simulation.restart_options - restarts_path = "simulation/restarts/" - restart_file_path = "phare_outputs" - - if "dir" in restart_options: - restart_file_path = restart_options["dir"] - - if "restart_time" in restart_options: - from pyphare.cpp import cpp_etc_lib - - restart_time = restart_options["restart_time"] - restart_file_load_path = cpp_etc_lib().restart_path_for_time( - restart_file_path, restart_time - ) - - if not os.path.exists(restart_file_load_path): - raise ValueError( - f"PHARE restart file not found for time {restart_time}" - ) - - deserialized_simulation = deserialize_sim( - _serialized_simulation_string(restart_file_load_path) - ) - if not simulation.is_restartable_compared_to(deserialized_simulation): - raise ValueError( - "deserialized Restart simulation is incompatible with configured simulation parameters" - ) - - add_vector_int( - restarts_path + "restart_ids", _patch_data_ids(restart_file_load_path) - ) - add_string(restarts_path + "loadPath", restart_file_load_path) - add_double(restarts_path + "restart_time", restart_time) - - if "mode" in restart_options: - add_string(restarts_path + "mode", restart_options["mode"]) - - add_string(restarts_path + "filePath", restart_file_path) - - if "elapsed_timestamps" in restart_options: - pp.add_array_as_vector( - restarts_path + "elapsed_timestamps", - restart_options["elapsed_timestamps"], - ) - - if "timestamps" in restart_options: - pp.add_array_as_vector( - restarts_path + "write_timestamps", restart_options["timestamps"] - ) + from .global_vars import sim + from . import initialize - add_string(restarts_path + "serialized_simulation", serialized_sim) - #### restarts added + initialize.general.populateDict(sim) - #### adding electrons - if simulation.electrons is None: - raise RuntimeError("Error - no electrons registered to this Simulation") + if sim.init_options is None: + initialize.user_fns.populateDict(sim) else: - for item in simulation.electrons.dict_path(): - if isinstance(item[1], str): - add_string("simulation/" + item[0], item[1]) - else: - add_double("simulation/" + item[0], item[1]) + initialize.samrai_hdf5.populateDict(sim) diff --git a/pyphare/pyphare/pharein/initialize/__init__.py b/pyphare/pyphare/pharein/initialize/__init__.py new file mode 100644 index 000000000..bf53e1c00 --- /dev/null +++ b/pyphare/pyphare/pharein/initialize/__init__.py @@ -0,0 +1,9 @@ +from . import general +from . import user_fns +from . import samrai_hdf5 + +__all__ = [ + "general", + "user_fns", + "samrai_hdf5", +] diff --git a/pyphare/pyphare/pharein/initialize/general.py b/pyphare/pyphare/pharein/initialize/general.py new file mode 100644 index 000000000..0bff869e2 --- /dev/null +++ b/pyphare/pyphare/pharein/initialize/general.py @@ -0,0 +1,318 @@ +# +# + +import os +import numpy as np +import pybindlibs.dictator as pp +from pyphare.core.phare_utilities import is_scalar + +from pyphare.pharein.load_balancer import LoadBalancer +from pyphare.pharein.simulation import ( + serialize as serialize_sim, + deserialize as deserialize_sim, +) + + +def _patch_data_ids(restart_file_dir): + """ + for restarts we save samrai patch data ids to the restart files, which we access from here + to tell samrai which patch datas to load from the restart file on restart + """ + from pyphare.cpp import cpp_etc_lib + + return cpp_etc_lib().patch_data_ids(restart_file_dir) + + +def _serialized_simulation_string(restart_file_dir): + from pyphare.cpp import cpp_etc_lib + + return cpp_etc_lib().serialized_simulation_string(restart_file_dir) + + +# converts scalars to array of expected size +# converts lists to arrays +class py_fn_wrapper: + def __init__(self, fn): + self.fn = fn + + def __call__(self, *xyz): + args = [np.asarray(arg) for arg in xyz] + ret = self.fn(*args) + if isinstance(ret, list): + ret = np.asarray(ret) + if is_scalar(ret): + ret = np.full(len(args[-1]), ret) + return ret + + +# Wrap calls to user init functions to turn C++ vectors to ndarrays, +# and returned ndarrays to C++ span +class fn_wrapper(py_fn_wrapper): + def __init__(self, fn): + super().__init__(fn) + + def __call__(self, *xyz): + from pyphare.cpp import cpp_etc_lib + + # convert numpy array to C++ SubSpan + # couples vector init functions to C++ + return cpp_etc_lib().makePyArrayWrapper(super().__call__(*xyz)) + + +# pybind complains if receiving wrong type +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)) + + +def add_size_t(path, val): + casted = int(val) + if casted < 0: + raise RuntimeError("pyphare.__init__::add_size_t received negative value") + pp.add_size_t(path, casted) + + +def add_vector_int(path, val): + pp.add_vector_int(path, list(val)) + + +add_string = pp.add_string + + +def populateDict(sim): + import pybindlibs.dictator as pp + + add_string("simulation/name", "simulation_test") + add_int("simulation/dimension", sim.ndim) + + if sim.smallest_patch_size is not None: + add_vector_int("simulation/AMR/smallest_patch_size", sim.smallest_patch_size) + if sim.largest_patch_size is not None: + add_vector_int("simulation/AMR/largest_patch_size", sim.largest_patch_size) + + add_string("simulation/grid/layout_type", sim.layout) + add_int("simulation/grid/nbr_cells/x", sim.cells[0]) + add_double("simulation/grid/meshsize/x", sim.dl[0]) + add_double("simulation/grid/origin/x", sim.origin[0]) + add_string("simulation/grid/boundary_type/x", sim.boundary_types[0]) + + if sim.ndim > 1: + add_int("simulation/grid/nbr_cells/y", sim.cells[1]) + add_double("simulation/grid/meshsize/y", sim.dl[1]) + add_double("simulation/grid/origin/y", sim.origin[1]) + add_string("simulation/grid/boundary_type/y", sim.boundary_types[1]) + + if sim.ndim > 2: + add_int("simulation/grid/nbr_cells/z", sim.cells[2]) + add_double("simulation/grid/meshsize/z", sim.dl[2]) + add_double("simulation/grid/origin/z", sim.origin[2]) + add_string("simulation/grid/boundary_type/z", sim.boundary_types[2]) + + add_int("simulation/interp_order", sim.interp_order) + add_int("simulation/refined_particle_nbr", sim.refined_particle_nbr) + add_double("simulation/time_step", sim.time_step) + add_int("simulation/time_step_nbr", sim.time_step_nbr) + + add_string("simulation/AMR/clustering", sim.clustering) + add_int("simulation/AMR/max_nbr_levels", sim.max_nbr_levels) + add_vector_int("simulation/AMR/nesting_buffer", sim.nesting_buffer) + + add_int("simulation/AMR/tag_buffer", sim.tag_buffer) + + refinement_boxes = sim.refinement_boxes + + def as_paths(rb): + add_int("simulation/AMR/refinement/boxes/nbr_levels/", len(rb.keys())) + for level, boxes in rb.items(): + level_path = "simulation/AMR/refinement/boxes/" + level + "/" + add_int(level_path + "nbr_boxes/", int(len(boxes))) + for box_i, box in enumerate(boxes): + box_id = "B" + str(box_i) + lower = box.lower + upper = box.upper + box_lower_path_x = box_id + "/lower/x/" + box_upper_path_x = box_id + "/upper/x/" + add_int(level_path + box_lower_path_x, lower[0]) + add_int(level_path + box_upper_path_x, upper[0]) + if len(lower) >= 2: + box_lower_path_y = box_id + "/lower/y/" + box_upper_path_y = box_id + "/upper/y/" + add_int(level_path + box_lower_path_y, lower[1]) + add_int(level_path + box_upper_path_y, upper[1]) + if len(lower) == 3: + box_lower_path_z = box_id + "/lower/z/" + box_upper_path_z = box_id + "/upper/z/" + add_int(level_path + box_lower_path_z, lower[2]) + add_int(level_path + box_upper_path_z, upper[2]) + + if refinement_boxes is not None and sim.refinement == "boxes": + as_paths(refinement_boxes) + elif sim.refinement == "tagging": + add_string("simulation/AMR/refinement/tagging/method", "auto") + # the two following params are hard-coded for now + # they will become configurable when we have multi-models or several methods + # per model + add_string("simulation/AMR/refinement/tagging/model", "HybridModel") + add_string("simulation/AMR/refinement/tagging/method", "default") + add_double("simulation/AMR/refinement/tagging/threshold", sim.tagging_threshold) + else: + add_string( + "simulation/AMR/refinement/tagging/method", "none" + ) # integrator.h might want some looking at + + add_string("simulation/algo/ion_updater/pusher/name", sim.particle_pusher) + add_double("simulation/algo/ohm/resistivity", sim.resistivity) + add_string("simulation/algo/ohm/hyper_mode", sim.hyper_mode) + add_double("simulation/algo/ohm/hyper_resistivity", sim.hyper_resistivity) + + # load balancer block start + lb = sim.load_balancer or LoadBalancer(active=False, _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 = sim.model + + if init_model.nbr_populations() < 0: + raise RuntimeError("Number of populations cannot be negative") + add_size_t("simulation/ions/nbrPopulations", init_model.nbr_populations()) + + modelDict = init_model.model_dict + for pop_index, pop in enumerate(init_model.populations): + pop_path = "simulation/ions/pop" + add_string(pop_path + "{:d}/name".format(pop_index), pop) + d = modelDict[pop] + add_double(pop_path + "{:d}/mass".format(pop_index), d["mass"]) + + add_string("simulation/electromag/name", "EM") + add_string("simulation/electromag/electric/name", "E") + add_string("simulation/electromag/magnetic/name", "B") + + serialized_sim = serialize_sim(sim) + + #### adding diagnostics + + diag_path = "simulation/diagnostics/" + for diag in list(sim.diagnostics.values()): + diag.attributes["serialized_simulation"] = serialized_sim + + type_path = diag_path + diag.type + "/" + name_path = type_path + diag.name + add_string(name_path + "/" + "type", diag.type) + add_string(name_path + "/" + "quantity", diag.quantity) + add_size_t(name_path + "/" + "flush_every", diag.flush_every) + pp.add_array_as_vector( + name_path + "/" + "write_timestamps", diag.write_timestamps + ) + pp.add_array_as_vector( + name_path + "/" + "compute_timestamps", diag.compute_timestamps + ) + + add_size_t(name_path + "/" + "n_attributes", len(diag.attributes)) + for attr_idx, attr_key in enumerate(diag.attributes): + add_string(name_path + "/" + f"attribute_{attr_idx}_key", attr_key) + add_string( + name_path + "/" + f"attribute_{attr_idx}_value", + diag.attributes[attr_key], + ) + + if len(sim.diagnostics) > 0: + if sim.diag_options is not None and "options" in sim.diag_options: + add_string(diag_path + "filePath", sim.diag_options["options"]["dir"]) + if "mode" in sim.diag_options["options"]: + add_string(diag_path + "mode", sim.diag_options["options"]["mode"]) + if "fine_dump_lvl_max" in sim.diag_options["options"]: + add_int( + diag_path + "fine_dump_lvl_max", + sim.diag_options["options"]["fine_dump_lvl_max"], + ) + else: + add_string(diag_path + "filePath", "phare_output") + #### diagnostics added + + #### adding restarts + if sim.restart_options is not None: + restart_options = sim.restart_options + restarts_path = "simulation/restarts/" + restart_file_path = "phare_outputs" + + if "dir" in restart_options: + restart_file_path = restart_options["dir"] + + if "restart_time" in restart_options: + from pyphare.cpp import cpp_etc_lib + + restart_time = restart_options["restart_time"] + restart_file_load_path = cpp_etc_lib().restart_path_for_time( + restart_file_path, restart_time + ) + + if not os.path.exists(restart_file_load_path): + raise ValueError( + f"PHARE restart file not found for time {restart_time}" + ) + + deserialized_simulation = deserialize_sim( + _serialized_simulation_string(restart_file_load_path) + ) + if not sim.is_restartable_compared_to(deserialized_simulation): + raise ValueError( + "deserialized Restart simulation is incompatible with configured simulation parameters" + ) + + add_vector_int( + restarts_path + "restart_ids", _patch_data_ids(restart_file_load_path) + ) + add_string(restarts_path + "loadPath", restart_file_load_path) + add_double(restarts_path + "restart_time", restart_time) + + if "mode" in restart_options: + add_string(restarts_path + "mode", restart_options["mode"]) + + add_string(restarts_path + "filePath", restart_file_path) + + if "elapsed_timestamps" in restart_options: + pp.add_array_as_vector( + restarts_path + "elapsed_timestamps", + restart_options["elapsed_timestamps"], + ) + + if "timestamps" in restart_options: + pp.add_array_as_vector( + restarts_path + "write_timestamps", restart_options["timestamps"] + ) + + add_string(restarts_path + "serialized_simulation", serialized_sim) + #### restarts added + + #### adding electrons + if sim.electrons is None: + raise RuntimeError("Error - no electrons registered to this Simulation") + else: + for item in sim.electrons.dict_path(): + if isinstance(item[1], str): + add_string("simulation/" + item[0], item[1]) + else: + add_double("simulation/" + item[0], item[1]) diff --git a/pyphare/pyphare/pharein/initialize/samrai_hdf5.py b/pyphare/pyphare/pharein/initialize/samrai_hdf5.py new file mode 100644 index 000000000..9c81afa5b --- /dev/null +++ b/pyphare/pyphare/pharein/initialize/samrai_hdf5.py @@ -0,0 +1,16 @@ +from .general import add_string, add_int + +# import pybindlibs.dictator as pp + + +def populateDict(sim): + init_model = sim.model + partinit = "particle_initializer" + for pop_index, pop in enumerate(init_model.populations): + pop_path = "simulation/ions/pop" + partinit_path = pop_path + "{:d}/".format(pop_index) + partinit + "/" + + add_string(partinit_path + "name", "samraih5") + add_string(partinit_path + "filepath", sim.init_options["dir"]) + add_int(partinit_path + "mpi_size", sim.init_options.get("mpi_size", 1)) + add_int(partinit_path + "index", sim.init_options.get("index", 0)) diff --git a/pyphare/pyphare/pharein/initialize/user_fns.py b/pyphare/pyphare/pharein/initialize/user_fns.py new file mode 100644 index 000000000..855f127e2 --- /dev/null +++ b/pyphare/pyphare/pharein/initialize/user_fns.py @@ -0,0 +1,51 @@ +from .general import add_double, add_string, add_int, add_size_t, fn_wrapper +import pybindlibs.dictator as pp + + +def populateDict(sim): + populate_electromag(sim) + populate_particles(sim) + + +def populate_electromag(sim): + addInitFunction = getattr(pp, "addInitFunction{:d}".format(sim.ndim) + "D") + modelDict = sim.model.model_dict + maginit_path = "simulation/electromag/magnetic/initializer/" + addInitFunction(maginit_path + "x_component", fn_wrapper(modelDict["bx"])) + addInitFunction(maginit_path + "y_component", fn_wrapper(modelDict["by"])) + addInitFunction(maginit_path + "z_component", fn_wrapper(modelDict["bz"])) + + +def populate_particles(sim): + addInitFunction = getattr(pp, "addInitFunction{:d}".format(sim.ndim) + "D") + + init_model = sim.model + modelDict = init_model.model_dict + + if init_model.nbr_populations() < 0: + raise RuntimeError("Number of populations cannot be negative") + add_size_t("simulation/ions/nbrPopulations", init_model.nbr_populations()) + + partinit = "particle_initializer" + for pop_index, pop in enumerate(init_model.populations): + pop_path = "simulation/ions/pop" + partinit_path = pop_path + "{:d}/".format(pop_index) + partinit + "/" + d = modelDict[pop] + + add_string(partinit_path + "name", "maxwellian") + + addInitFunction(partinit_path + "density", fn_wrapper(d["density"])) + addInitFunction(partinit_path + "bulk_velocity_x", fn_wrapper(d["vx"])) + addInitFunction(partinit_path + "bulk_velocity_y", fn_wrapper(d["vy"])) + addInitFunction(partinit_path + "bulk_velocity_z", fn_wrapper(d["vz"])) + 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_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"]) diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 5fb41ecc3..17bc2ac1c 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -534,6 +534,23 @@ def check_restart_options(**kwargs): return restart_options +def check_init_options(**kwargs): + """Advanced options to initialize from SAMRAI HDF5 files""" + + formats = ["samraih5"] + init_options = kwargs.get("init_options", None) + + if init_options is not None and "format" in init_options: + if init_options["format"] not in formats: + raise ValueError("Error - init_options format is invalid") + if "options" in init_options and "dir" in init_options["options"]: + init_options["options"]["dir"] = check_directory( + init_options["options"]["dir"], "init_options" + ) + + return init_options + + def validate_restart_options(sim): import pyphare.pharein.restarts as restarts @@ -658,6 +675,7 @@ def wrapper(simulation_object, **kwargs): "description", "dry_run", "write_reports", + "init_options", ] accepted_keywords += check_optional_keywords(**kwargs) @@ -692,6 +710,7 @@ def wrapper(simulation_object, **kwargs): ndim = compute_dimension(cells) kwargs["diag_options"] = check_diag_options(**kwargs) + kwargs["init_options"] = check_init_options(**kwargs) kwargs["restart_options"] = check_restart_options(**kwargs) kwargs["boundary_types"] = check_boundaries(ndim, **kwargs) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index 487dc8605..eb47ce4d5 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -1,6 +1,7 @@ -from dataclasses import dataclass, field -from copy import deepcopy +import os import numpy as np +from copy import deepcopy +from dataclasses import dataclass, field from typing import Any, List, Tuple @@ -598,14 +599,16 @@ def __bool__(self): return not self.failed def __repr__(self): + strs = [] for msg, ref, cmp in self: - print(msg) + strs.append(msg) try: if type(ref) is FieldData: phut.assert_fp_any_all_close(ref[:], cmp[:], atol=1e-16) except AssertionError as e: - print(e) - return self.failed[0][0] + strs.append(e) + strs.append(self.failed[0][0]) + return os.linesep.join(strs) def __call__(self, reason, ref=None, cmp=None): self.failed.append((reason, ref, cmp)) @@ -705,6 +708,7 @@ def _skip(qty): l0_pds[k][patch.box] = v[patch.box] elif isinstance(v, ParticleData): l0_pds[k].dataset.add(v.dataset) + else: raise RuntimeError("unexpected state") return cier diff --git a/pyphare/pyphare/pharesee/hierarchy/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/patchdata.py index e0727b9a2..b8b9ba815 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patchdata.py +++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py @@ -1,6 +1,8 @@ import numpy as np + from ...core import phare_utilities as phut + from ...core import box as boxm from ...core.box import Box diff --git a/res/cmake/test.cmake b/res/cmake/test.cmake index 1f7331b1d..c3b81bd1d 100644 --- a/res/cmake/test.cmake +++ b/res/cmake/test.cmake @@ -16,7 +16,6 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests add_subdirectory(tests/core/data/electrons) add_subdirectory(tests/core/data/ion_population) add_subdirectory(tests/core/data/maxwellian_particle_initializer) - add_subdirectory(tests/core/data/particle_initializer) add_subdirectory(tests/core/utilities/box) add_subdirectory(tests/core/utilities/range) add_subdirectory(tests/core/utilities/index) @@ -42,6 +41,7 @@ if (test AND ${PHARE_EXEC_LEVEL_MIN} GREATER 0) # 0 = no tests add_subdirectory(tests/amr/data/field/refine) add_subdirectory(tests/amr/data/field/variable) add_subdirectory(tests/amr/data/field/time_interpolate) + add_subdirectory(tests/amr/data/particles/initializer) add_subdirectory(tests/amr/resources_manager) add_subdirectory(tests/amr/messengers) add_subdirectory(tests/amr/models) diff --git a/src/amr/data/electromag/electromag_initializer.hpp b/src/amr/data/electromag/electromag_initializer.hpp new file mode 100644 index 000000000..3dab99726 --- /dev/null +++ b/src/amr/data/electromag/electromag_initializer.hpp @@ -0,0 +1,92 @@ +#ifndef _PHARE_AMR_DATA_ELECTROMAG_ELECTROMAG_INITIALIZER_HPP_ +#define _PHARE_AMR_DATA_ELECTROMAG_ELECTROMAG_INITIALIZER_HPP_ + +#include "core/data/vecfield/vecfield_initializer.hpp" +#include "amr/data/field/initializers/samrai_hdf5_field_initializer.hpp" +#include "initializer/data_provider.hpp" + + +#include + +namespace PHARE::amr +{ + +template +class ElectromagInitializer +{ +public: + virtual void init(Electromag_t& /*em*/, GridLayout const& /*layout*/) const { /*noop*/ } + virtual ~ElectromagInitializer() {} +}; + +template +class ElectromagUserFuncInitializer : public ElectromagInitializer +{ + using VecFieldInit = core::VecFieldInitializer; + +public: + ElectromagUserFuncInitializer(initializer::PHAREDict const& dict) + : Binit_{dict["magnetic"]["initializer"]} + { + } + virtual ~ElectromagUserFuncInitializer() {} + + void init(Electromag_t& em, GridLayout const& layout) const override + { + Binit_.initialize(em.B, layout); + } + + VecFieldInit Binit_; +}; + + +/* +/PHARE_hierarchy/level_0000/level_0000-patch_0000000-block_0000000/EM_B_y##default/d_box +/PHARE_hierarchy/level_0000/level_0000-patch_0000000-block_0000000/EM_B_x##default/field_EM_B_x +/PHARE_hierarchy/level_0000/level_0000-patch_0000000-block_0000000/EM_B_y##default/field_EM_B_y +/PHARE_hierarchy/level_0000/level_0000-patch_0000000-block_0000000/EM_B_z##default/field_EM_B_z +*/ + + +template +class ElectromagSamraiH5Initializer : public ElectromagInitializer +{ + using vecfield_type = typename Electromag_t::vecfield_type; + using field_type = typename vecfield_type::field_type; + + +public: + ElectromagSamraiH5Initializer(initializer::PHAREDict const& dict) + : dict_{dict} + { + } + virtual ~ElectromagSamraiH5Initializer() {} + + void init(Electromag_t& em, GridLayout const& layout) const override + { + for (auto& field : em.B) + SamraiHDF5FieldInitializer{}.load(field, layout); + } + + initializer::PHAREDict const dict_; +}; + + + +class ElectromagInitializerFactory +{ +public: + template + NO_DISCARD static std::unique_ptr> + create(initializer::PHAREDict const& dict) + { + if (dict["magnetic"]["initializer"].contains("x_component")) + return std::make_unique>(dict); + else + return std::make_unique>(dict); + } +}; + +} // namespace PHARE::amr + +#endif // _PHARE_AMR_DATA_ELECTROMAG_ELECTROMAG_INITIALIZER_HPP_ diff --git a/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp b/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp new file mode 100644 index 000000000..ffce2ab5c --- /dev/null +++ b/src/amr/data/field/initializers/samrai_hdf5_field_initializer.hpp @@ -0,0 +1,69 @@ +#ifndef _PHARE_AMR_DATA_FIELD_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ +#define _PHARE_AMR_DATA_FIELD_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ + +#include + +#include "core/data/grid/gridlayoutdefs.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" + +#include "core/utilities/types.hpp" +#include "core/utilities/point/point.hpp" + +#include "amr/data/initializers/samrai_hdf5_initializer.hpp" + +namespace PHARE::amr +{ + + +template +class SamraiHDF5FieldInitializer // no virtual classes needed (yet) +{ +public: + static constexpr auto dimension = GridLayout::dimension; + + SamraiHDF5FieldInitializer() {} + + void load(Field_t& Field, GridLayout const& layout) const; +}; + + +template +void SamraiHDF5FieldInitializer::load(Field_t& field, + GridLayout const& layout) const +{ + auto const& dest_box = layout.AMRBox(); + auto const& centering = layout.centering(field.physicalQuantity()); + auto const& overlaps = SamraiH5Interface::INSTANCE().box_intersections(dest_box); + for (auto const& [overlap_box, h5FilePtr, pdataptr] : overlaps) + { + auto& h5File = *h5FilePtr; + auto& pdata = *pdataptr; + auto const src_box = pdata.box; + auto const data = h5File.template read_data_set_flat( + pdata.base_path + "/" + field.name() + "##default/field_" + field.name()); + core::Box const lcl_src_gbox{ + core::Point{core::ConstArray()}, + core::Point{ + core::for_N([&](auto i) { + return static_cast( + src_box.upper[i] - src_box.lower[i] + (GridLayout::nbrGhosts() * 2) + + (centering[i] == core::QtyCentering::primal ? 1 : 0)); + })}}; + auto const data_view = core::make_array_view(data.data(), *lcl_src_gbox.shape()); + auto const overlap_gb = grow(overlap_box, GridLayout::nbrGhosts()); + auto const lcl_src_box = layout.AMRToLocal(overlap_gb, src_box); + auto const lcl_dst_box = layout.AMRToLocal(overlap_gb, dest_box); + auto src_it = lcl_src_box.begin(); + auto dst_it = lcl_dst_box.begin(); + for (; src_it != lcl_src_box.end(); ++src_it, ++dst_it) + field(*dst_it) = data_view(*src_it); + } +} + + + +} // namespace PHARE::amr + + + +#endif diff --git a/src/amr/data/initializers/samrai_hdf5_initializer.hpp b/src/amr/data/initializers/samrai_hdf5_initializer.hpp new file mode 100644 index 000000000..3b7def9d1 --- /dev/null +++ b/src/amr/data/initializers/samrai_hdf5_initializer.hpp @@ -0,0 +1,140 @@ +#ifndef _PHARE_AMR_DATA_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ +#define _PHARE_AMR_DATA_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ + +#include +#include + +#include "core/def.hpp" +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "hdf5/detail/h5/h5_file.hpp" + +#include "SAMRAI/tbox/HDFDatabase.h" + + +namespace PHARE::amr +{ + +template +struct SamraiH5PatchDataInfo +{ + using Box_t = core::Box; + + SamraiH5PatchDataInfo(Box_t const& b, std::string const& p) + : box{b} + , base_path{p} + { + } + + + Box_t const box; + std::string const base_path; +}; + + + +template +class SamraiH5Interface +{ + struct SamraiHDF5File; + +public: + using Box_t = core::Box; + + static SamraiH5Interface& INSTANCE() + { + static SamraiH5Interface i; + return i; + } + + void populate_from(std::string const& dir, int const& idx, int const& mpi_size, + std::string const& field_name = "field_EM_B_x"); + + NO_DISCARD auto static getRestartFileFullPath(std::string path, int const& idx, + int const& mpi_size, int const& rank) + { + return path // + + "/restore." + SAMRAI::tbox::Utilities::intToString(idx, 6) // + + "/nodes." + SAMRAI::tbox::Utilities::nodeToString(mpi_size) // + + "/proc." + SAMRAI::tbox::Utilities::processorToString(rank); + } + + auto box_intersections(Box_t const& box) + { + using Tup = std::tuple const* const>; + std::vector overlaps; + for (auto const& h5File : restart_files) + for (auto const& patch : h5File->patches) + if (auto const intersection = box * patch.box) + overlaps.emplace_back(*intersection, h5File.get(), &patch); + return overlaps; + } + + + +private: + std::vector> restart_files; + std::unordered_map box2dataset; +}; + + + +template +struct SamraiH5Interface::SamraiHDF5File : public hdf5::h5::HighFiveFile +{ + using Super = hdf5::h5::HighFiveFile; + + SamraiHDF5File(std::string const& filepath) + : Super{filepath, HighFive::File::ReadOnly, /*para=*/false} + , filepath_{filepath} + { + } + + auto getBoxFromPath(std::string const& path) const + { + std::size_t constexpr samrai_dim = 3; // always 3! + auto constexpr _to_std_array = [](auto& i) { + return core::sized_array( + *reinterpret_cast const*>(&i)); + }; + + SAMRAI::tbox::HDFDatabase db{"db"}; + db.open(filepath_); + auto const boxes = db.getDatabaseBoxVector(path); + return Box_t{_to_std_array(boxes[0].d_data.d_lo), _to_std_array(boxes[0].d_data.d_hi)}; + } + + std::string const filepath_; + std::vector> patches; +}; + + + +template +void SamraiH5Interface::populate_from(std::string const& dir, int const& idx, + int const& mpi_size, + std::string const& field_name) +{ + if (restart_files.size()) // executed per pop, but we only need to run this once + return; + + for (int rank = 0; rank < mpi_size; ++rank) + { + auto const hdf5_filepath = getRestartFileFullPath(dir, idx, mpi_size, rank); + auto& h5File = *restart_files.emplace_back(std::make_unique(hdf5_filepath)); + for (auto const& group : h5File.scan_for_groups({"level_0000", field_name})) + { + auto const field_path = group.substr(0, group.rfind("/")); + auto const& field_box = h5File.getBoxFromPath(field_path + "/d_box"); + h5File.patches.emplace_back(field_box, field_path.substr(0, field_path.rfind("/"))); + } + } +} + + + +} // namespace PHARE::amr + + +#endif diff --git a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp b/src/amr/data/particles/initializers/particle_initializer_factory.hpp similarity index 70% rename from src/core/data/ions/particle_initializers/particle_initializer_factory.hpp rename to src/amr/data/particles/initializers/particle_initializer_factory.hpp index 97cfa6b0c..48681fdb7 100644 --- a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp +++ b/src/amr/data/particles/initializers/particle_initializer_factory.hpp @@ -5,19 +5,23 @@ #include "core/def.hpp" #include "core/utilities/types.hpp" #include "initializer/data_provider.hpp" -#include "maxwellian_particle_initializer.hpp" -#include "particle_initializer.hpp" + +#include "core/data/ions/particle_initializers/particle_initializer.hpp" +#include "core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp" + +#include "samrai_hdf5_particle_initializer.hpp" +#include "amr/data/initializers/samrai_hdf5_initializer.hpp" #include namespace PHARE { -namespace core +namespace amr { template class ParticleInitializerFactory { - using ParticleInitializerT = ParticleInitializer; + using ParticleInitializerT = core::ParticleInitializer; static constexpr auto dimension = GridLayout::dimension; @@ -66,8 +70,8 @@ namespace core if (basisName == "cartesian") { return std::make_unique< - MaxwellianParticleInitializer>( - density, v, vth, charge, nbrPartPerCell, seed, Basis::Cartesian, + core::MaxwellianParticleInitializer>( + density, v, vth, charge, nbrPartPerCell, seed, core::Basis::Cartesian, magneticField, densityCutOff); } else if (basisName == "magnetic") @@ -77,17 +81,30 @@ namespace core magneticField[2] = dict["magnetic_x"].template to(); return std::make_unique< - MaxwellianParticleInitializer>( - density, v, vth, charge, nbrPartPerCell, seed, Basis::Magnetic, + core::MaxwellianParticleInitializer>( + density, v, vth, charge, nbrPartPerCell, seed, core::Basis::Magnetic, magneticField, densityCutOff); } } - // TODO throw? - return nullptr; + + if (initializerName == "samraih5") + { + auto const dir = dict["filepath"].template to(); + int const index = dict["index"].template to(); + int const mpi_size = dict["mpi_size"].template to(); + + // scan restart files for later use + SamraiH5Interface::INSTANCE().populate_from(dir, index, mpi_size); + + return std::make_unique>(); + } + + + throw std::runtime_error("No Particle Initializer chosen!"); } }; -} // namespace core +} // namespace amr } // namespace PHARE diff --git a/src/amr/data/particles/initializers/samrai_hdf5_particle_initializer.hpp b/src/amr/data/particles/initializers/samrai_hdf5_particle_initializer.hpp new file mode 100644 index 000000000..197bad3d2 --- /dev/null +++ b/src/amr/data/particles/initializers/samrai_hdf5_particle_initializer.hpp @@ -0,0 +1,72 @@ +#ifndef _PHARE_AMR_DATA_PARTICLE_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ +#define _PHARE_AMR_DATA_PARTICLE_INITIAZILIZERS_SAMRAI_HDF5_INITIALIZER_HPP_ + +#include + +#include "core/utilities/types.hpp" +#include "core/utilities/box/box.hpp" +#include "core/utilities/point/point.hpp" +#include "core/data/ions/particle_initializers/particle_initializer.hpp" + +#include "core/data/particles/particle_packer.hpp" +#include "amr/data/initializers/samrai_hdf5_initializer.hpp" + + +namespace PHARE::amr +{ + + +template +class SamraiHDF5ParticleInitializer : public core::ParticleInitializer +{ +public: + static constexpr auto dimension = GridLayout::dimension; + + + + SamraiHDF5ParticleInitializer() {} + + + + void loadParticles(ParticleArray& particles, GridLayout const& layout, + std::string const& popname) const override; +}; + + + +template +void SamraiHDF5ParticleInitializer::loadParticles( + ParticleArray& particles, GridLayout const& layout, std::string const& popname) const +{ + using Packer = core::ParticlePacker; + + auto const& dest_box = layout.AMRBox(); + auto const& overlaps = SamraiH5Interface::INSTANCE().box_intersections(dest_box); + + for (auto const& [overlap_box, h5FilePtr, pdataptr] : overlaps) + { + auto& h5File = *h5FilePtr; + auto& pdata = *pdataptr; + + std::string const poppath = pdata.base_path + "/" + popname + "##default/domainParticles_"; + core::ContiguousParticles soa{0}; + + { + std::size_t part_idx = 0; + core::apply(soa.as_tuple(), [&](auto& arg) { + auto const datapath = poppath + Packer::keys()[part_idx++]; + h5File.file().getDataSet(datapath).read(arg); + }); + } + + for (std::size_t i = 0; i < soa.size(); ++i) + if (auto const p = soa.copy(i); core::isIn(core::Point{p.iCell}, dest_box)) + particles.push_back(p); + } +} + + +} // namespace PHARE::amr + + +#endif diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index e0ef8386a..3acdf4fcd 100644 --- a/src/amr/level_initializer/hybrid_level_initializer.hpp +++ b/src/amr/level_initializer/hybrid_level_initializer.hpp @@ -4,11 +4,9 @@ #include "amr/level_initializer/level_initializer.hpp" #include "amr/messengers/hybrid_messenger.hpp" #include "amr/messengers/messenger.hpp" -#include "amr/physical_models/hybrid_model.hpp" #include "amr/physical_models/physical_model.hpp" #include "amr/resources_manager/amr_utils.hpp" #include "core/data/grid/gridlayout_utils.hpp" -#include "core/data/ions/ions.hpp" #include "core/numerics/ampere/ampere.hpp" #include "core/numerics/interpolator/interpolator.hpp" #include "core/numerics/moments/moments.hpp" @@ -43,10 +41,12 @@ namespace solver : ohm_{dict["algo"]["ohm"]} { } - virtual void initialize(std::shared_ptr const& hierarchy, int levelNumber, - std::shared_ptr const& oldLevel, IPhysicalModelT& model, - amr::IMessenger& messenger, double initDataTime, - bool isRegridding) override + + + void initialize(std::shared_ptr const& hierarchy, int levelNumber, + std::shared_ptr const& oldLevel, IPhysicalModelT& model, + amr::IMessenger& messenger, double initDataTime, + bool isRegridding) override { core::Interpolator interpolate_; auto& hybridModel = static_cast(model); @@ -163,6 +163,8 @@ namespace solver hybMessenger.prepareStep(hybridModel, level, initDataTime); } }; + + } // namespace solver } // namespace PHARE diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index 449bafc22..7488ab021 100644 --- a/src/amr/physical_models/hybrid_model.hpp +++ b/src/amr/physical_models/hybrid_model.hpp @@ -1,16 +1,23 @@ #ifndef PHARE_HYBRID_MODEL_HPP #define PHARE_HYBRID_MODEL_HPP +#include #include -#include "initializer/data_provider.hpp" +#include "core/def.hpp" +#include "core/data/vecfield/vecfield.hpp" #include "core/models/hybrid_state.hpp" + +#include "initializer/data_provider.hpp" + #include "amr/physical_models/physical_model.hpp" -#include "core/data/ions/particle_initializers/particle_initializer_factory.hpp" +#include "amr/data/particles/initializers/particle_initializer_factory.hpp" + +#include "amr/data/electromag/electromag_initializer.hpp" + #include "amr/resources_manager/resources_manager.hpp" #include "amr/messengers/hybrid_messenger_info.hpp" -#include "core/data/vecfield/vecfield.hpp" -#include "core/def.hpp" + namespace PHARE::solver { @@ -22,6 +29,9 @@ template class HybridModel : public IPhysicalModel { + struct _Initializers; + + public: static constexpr auto dimension = GridLayoutT::dimension; @@ -41,12 +51,12 @@ class HybridModel : public IPhysicalModel using particle_array_type = typename Ions::particle_array_type; using resources_manager_type = amr::ResourcesManager; using ParticleInitializerFactory - = core::ParticleInitializerFactory; - + = amr::ParticleInitializerFactory; + using HybridState_t = core::HybridState; static const inline std::string model_name = "HybridModel"; - core::HybridState state; + HybridState_t state; std::shared_ptr resourcesManager; @@ -84,6 +94,7 @@ class HybridModel : public IPhysicalModel : IPhysicalModel{model_name} , state{dict} , resourcesManager{std::move(_resourcesManager)} + , initializers_{std::make_unique<_Initializers>(dict, state)} { } @@ -107,6 +118,9 @@ class HybridModel : public IPhysicalModel //------------------------------------------------------------------------- std::unordered_map>> tags; + std::unique_ptr<_Initializers> initializers_; + + auto& initializers() { return *initializers_; } }; @@ -130,17 +144,12 @@ void HybridModel::i auto _ = this->resourcesManager->setOnPatch(*patch, state.electromag, state.ions); for (auto& pop : ions) - { - auto const& info = pop.particleInitializerInfo(); - auto particleInitializer = ParticleInitializerFactory::create(info); - particleInitializer->loadParticles(pop.domainParticles(), layout); - } - - state.electromag.initialize(layout); + initializers_->particles(pop).loadParticles(pop.domainParticles(), layout, pop.name()); + initializers_->electromag().init(state.electromag, layout); } - resourcesManager->registerForRestarts(*this); + initializers_.release(); } @@ -191,6 +200,40 @@ struct type_list_to_hybrid_model template using type_list_to_hybrid_model_t = typename type_list_to_hybrid_model::type; + + +template +struct HybridModel::_Initializers +{ + using particle_array_type = typename Ions::particle_array_type; + using ParticleInitializer_t = core::ParticleInitializer; + using ParticleInitializerFactory + = amr::ParticleInitializerFactory; + using ElectromagInitializer_t = amr::ElectromagInitializer; + + _Initializers(initializer::PHAREDict const& dict, HybridState_t const& state) + : dict_{dict} + , electromag_init{amr::ElectromagInitializerFactory::create( + dict["electromag"])} + { + for (auto& pop : state.ions) + particle_pop_init.emplace( + pop.name(), ParticleInitializerFactory::create(pop.particleInitializerInfo())); + } + + auto& electromag() { return *electromag_init; } + + auto& particles(typename Ions::value_type& pop) { return *particle_pop_init[pop.name()]; } + + initializer::PHAREDict dict_; + std::unique_ptr electromag_init; + std::unordered_map> particle_pop_init; +}; + + + + } // namespace PHARE::solver #endif diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 293de3711..d9754f624 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -19,7 +19,6 @@ set( SOURCES_INC data/electrons/electrons.hpp data/ions/particle_initializers/particle_initializer.hpp data/ions/particle_initializers/maxwellian_particle_initializer.hpp - data/ions/particle_initializers/particle_initializer_factory.hpp data/tensorfield/tensorfield.hpp data/vecfield/vecfield.hpp data/vecfield/vecfield_component.hpp diff --git a/src/core/data/electromag/electromag.hpp b/src/core/data/electromag/electromag.hpp index ed0a766f3..880c12e00 100644 --- a/src/core/data/electromag/electromag.hpp +++ b/src/core/data/electromag/electromag.hpp @@ -4,10 +4,10 @@ #include #include +#include "core/def.hpp" #include "core/hybrid/hybrid_quantities.hpp" -#include "core/data/vecfield/vecfield_initializer.hpp" + #include "initializer/data_provider.hpp" -#include "core/def.hpp" namespace PHARE @@ -17,16 +17,19 @@ namespace core template class Electromag { + using This = Electromag; + public: - static constexpr std::size_t dimension = VecFieldT::dimension; + using vecfield_type = VecFieldT; + auto static constexpr dimension = VecFieldT::dimension; explicit Electromag(std::string name) : E{name + "_E", HybridQuantity::Vector::E} , B{name + "_B", HybridQuantity::Vector::B} - , Binit_{} { } + explicit Electromag(initializer::PHAREDict const& dict) : E{dict["name"].template to() + "_" + dict["electric"]["name"].template to(), @@ -34,18 +37,9 @@ namespace core , B{dict["name"].template to() + "_" + dict["magnetic"]["name"].template to(), HybridQuantity::Vector::B} - , Binit_{dict["magnetic"]["initializer"]} { } - using vecfield_type = VecFieldT; - - - template - void initialize(GridLayout const& layout) - { - Binit_.initialize(B, layout); - } //------------------------------------------------------------------------- @@ -76,10 +70,13 @@ namespace core VecFieldT E; VecFieldT B; - - private: - VecFieldInitializer Binit_; }; + + } // namespace core } // namespace PHARE + + + + #endif diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 74297a864..f06f008a0 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -803,7 +803,8 @@ namespace core * This method only deals with **cell** indexes. */ template - NO_DISCARD auto AMRToLocal(Point AMRPoint) const + NO_DISCARD auto AMRToLocal(Point const& AMRPoint, + Box const& localbox) const { static_assert(std::is_integral_v, "Error, must be MeshIndex (integral Point)"); Point localPoint; @@ -811,33 +812,40 @@ namespace core // any direction, it's the same because we want cells auto localStart = physicalStartIndex(QtyCentering::dual, Direction::X); - // for (auto i = 0u; i < dimension; ++i) { - int local = AMRPoint[i] - (AMRBox_.lower[i] - localStart); + int local = AMRPoint[i] - (localbox.lower[i] - localStart); assert(local >= 0); localPoint[i] = local; } return localPoint; } + template + NO_DISCARD auto AMRToLocal(Point const& AMRPoint) const + { + return AMRToLocal(AMRPoint, AMRBox_); + } + /** * @brief AMRToLocal returns the local Box associated with the given AMR one. * This method only deals with **cell** indexes. */ template - NO_DISCARD auto AMRToLocal(Box AMRBox) const + NO_DISCARD auto AMRToLocal(Box const& AMRBox, + Box const& localbox) const { static_assert(std::is_integral_v, "Error, must be MeshIndex (integral Point)"); - auto localBox = Box{}; - - localBox.lower = AMRToLocal(AMRBox.lower); - localBox.upper = AMRToLocal(AMRBox.upper); - - return localBox; + return Box{AMRToLocal(AMRBox.lower, localbox), + AMRToLocal(AMRBox.upper, localbox)}; } + template + NO_DISCARD auto AMRToLocal(Box const& AMRBox) const + { + return AMRToLocal(AMRBox, AMRBox_); + } template @@ -1171,8 +1179,28 @@ namespace core evalOnBox_(field, fn, indices); } + auto levelNumber() const { return levelNumber_; } + + template + auto domainBoxFor(Field const& field) const + { + return _BoxFor(field, [&](auto const& centering, auto const direction) { + return this->physicalStartToEnd(centering, direction); + }); + } + + template + auto ghostBoxFor(Field const& field) const + { + return _BoxFor(field, [&](auto const& centering, auto const direction) { + return this->ghostStartToEnd(centering, direction); + }); + } + + + private: template static void evalOnBox_(Field& field, Fn& fn, IndicesFn& startToEnd) @@ -1207,6 +1235,20 @@ namespace core } + template + auto _BoxFor(Field const& field, Fn startToEnd) const + { + constexpr auto directions = std::array{Direction::X, Direction::Y, Direction::Z}; + std::array lower, upper; + core::for_N([&](auto i) { + auto const [i0, i1] = startToEnd(field, directions[i]); + lower[i] = i0; + upper[i] = i1; + }); + return Box{lower, upper}; + } + + template auto StartToEndIndices_(Centering const& centering, StartToEnd const&& startToEnd, bool const includeEnd = false) const diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 2701ebdac..553f5c59d 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -15,7 +15,7 @@ #include "core/hybrid/hybrid_quantities.hpp" #include "core/data/vecfield/vecfield_component.hpp" #include "initializer/data_provider.hpp" -#include "particle_initializers/particle_initializer_factory.hpp" + #include "core/utilities/algorithm.hpp" namespace PHARE 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 49098af06..7ec7b174b 100644 --- a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp +++ b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp @@ -62,7 +62,8 @@ class MaxwellianParticleInitializer : public ParticleInitializer void MaxwellianParticleInitializer::loadParticles( - ParticleArray& particles, GridLayout const& layout) const + ParticleArray& particles, GridLayout const& layout, std::string const& /*popname*/) const { auto point = [](std::size_t i, auto const& indices) -> core::Point { if constexpr (dimension == 1) diff --git a/src/core/data/ions/particle_initializers/particle_initializer.hpp b/src/core/data/ions/particle_initializers/particle_initializer.hpp index d8a01a071..1326955b5 100644 --- a/src/core/data/ions/particle_initializers/particle_initializer.hpp +++ b/src/core/data/ions/particle_initializers/particle_initializer.hpp @@ -10,7 +10,9 @@ namespace core class ParticleInitializer { public: - virtual void loadParticles(ParticleArray& particles, GridLayout const& layout) const = 0; + virtual void loadParticles(ParticleArray& particles, GridLayout const& layout, + std::string const& popname) const + = 0; virtual ~ParticleInitializer() = default; }; diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index 7d52d1d01..72a521fe9 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -2,7 +2,9 @@ #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP #include "core/def.hpp" +#include "core/logger.hpp" #include +#include #include #include #include @@ -19,20 +21,29 @@ namespace PHARE::core template struct NdArrayViewer { + template typename Indexes, typename Index> + NO_DISCARD static std::size_t idx(NCells const& nCells, Indexes const& indexes) + + { + if constexpr (dim == 1) + return idx(nCells, indexes[0]); + + else if constexpr (dim == 2) + return idx(nCells, indexes[0], indexes[1]); + + else if constexpr (dim == 3) + return idx(nCells, indexes[0], indexes[1], indexes[2]); + } + template - NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, - Indexes const&... indexes) + NO_DISCARD static std::size_t idx(NCells const& nCells, Indexes const&... indexes) { auto params = std::forward_as_tuple(indexes...); static_assert(sizeof...(Indexes) == dim); - // static_assert((... && std::is_unsigned_v)); TODO : manage later if - // this test should be included if constexpr (dim == 1) { - auto i = std::get<0>(params); - - return data[i]; + return std::get<0>(params); } if constexpr (dim == 2) @@ -41,9 +52,9 @@ struct NdArrayViewer auto j = std::get<1>(params); if constexpr (c_ordering) - return data[j + i * nCells[1]]; + return j + i * nCells[1]; else - return data[i + j * nCells[0]]; + return i + j * nCells[0]; } if constexpr (dim == 3) @@ -53,25 +64,28 @@ struct NdArrayViewer auto k = std::get<2>(params); if constexpr (c_ordering) - return data[k + j * nCells[2] + i * nCells[1] * nCells[2]]; + return k + j * nCells[2] + i * nCells[1] * nCells[2]; else - return data[i + j * nCells[0] + k * nCells[1] * nCells[0]]; + return i + j * nCells[0] + k * nCells[1] * nCells[0]; } } + + template + NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, + Indexes const&... indexes) + { + return data[idx(nCells, indexes...)]; + } + template typename Indexes, typename Index> NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, Indexes const& indexes) { - if constexpr (dim == 1) - return at(data, nCells, indexes[0]); - - else if constexpr (dim == 2) - return at(data, nCells, indexes[0], indexes[1]); - - else if constexpr (dim == 3) - return at(data, nCells, indexes[0], indexes[1], indexes[2]); + auto const i = idx(nCells, indexes); + assert(i < product(nCells, std::size_t{1})); // debug bounds check + return data[i]; } }; diff --git a/src/core/data/vecfield/vecfield_initializer.hpp b/src/core/data/vecfield/vecfield_initializer.hpp index edeb83bc6..38584d8c1 100644 --- a/src/core/data/vecfield/vecfield_initializer.hpp +++ b/src/core/data/vecfield/vecfield_initializer.hpp @@ -27,7 +27,7 @@ namespace core template - void initialize(VecField& v, GridLayout const& layout) + void initialize(VecField& v, GridLayout const& layout) const { static_assert(GridLayout::dimension == VecField::dimension, "dimension mismatch between vecfield and gridlayout"); diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index 0ff5d1a8a..c070e59b7 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -181,6 +181,11 @@ class box_iterator return *this; } + bool operator==(box_iterator const& other) const + { + return box_ == other.box_ and index_ == other.index_; + } + bool operator!=(box_iterator const& other) const { return box_ != other.box_ or index_ != other.index_; diff --git a/src/core/utilities/types.hpp b/src/core/utilities/types.hpp index ee90667ad..6a9618d36 100644 --- a/src/core/utilities/types.hpp +++ b/src/core/utilities/types.hpp @@ -149,6 +149,15 @@ namespace core return to; } + template + NO_DISCARD std::array to_array(std::vector const& from) + { + std::array to{}; + for (std::size_t i = 0; i < to.size(); i++) + to[i] = from[i]; + return to; + } + template NO_DISCARD constexpr auto as_sized_array(Args&&... args) diff --git a/src/hdf5/detail/h5/group_scanner.hpp b/src/hdf5/detail/h5/group_scanner.hpp new file mode 100644 index 000000000..f62a3daee --- /dev/null +++ b/src/hdf5/detail/h5/group_scanner.hpp @@ -0,0 +1,73 @@ +#ifndef _PHARE_HDF5_DETAIL_H5_GROUP_SCANNER_HPP_ +#define _PHARE_HDF5_DETAIL_H5_GROUP_SCANNER_HPP_ + +#include "core/def.hpp" +#include "core/def/phare_mpi.hpp" + +#include "highfive/H5File.hpp" +#include "highfive/H5Easy.hpp" + +#include "core/utilities/types.hpp" +#include "core/utilities/mpi_utils.hpp" +#include "core/utilities/meta/meta_utilities.hpp" + +#include + +namespace PHARE::hdf5::h5 +{ + +template +struct GroupScanner +{ + GroupScanner(HighFiveFile_t const& h5_, std::vector const& contains_ = {}) + : h5{h5_} + , contains{contains_} + { + } + + + auto& scan(std::string const& from = "/") + { + // prevent double / at start of full path + scan(h5.file().getGroup(from), from == "/" ? "" : from); + return groups; + } + void scan(HighFive::Group const& group, std::string const& path) + { + for (auto const& node : group.listObjectNames()) + { + if (group.getObjectType(node) == HighFive::ObjectType::Group) + scan(group.getGroup(node), path + "/" + node); + else + { + auto fpath = path + "/" + node; + if (contains.size() == 0) + groups.insert(fpath); + else + { + bool cont = false; + for (auto const& c : contains) + if (fpath.find(c) == std::string::npos) + { + cont = true; + break; + } + if (cont) // next node in listObjectNames + continue; + groups.insert(fpath); + } + } + } + } + + + HighFiveFile_t const& h5; + std::vector contains; + std::unordered_set groups{}; +}; + + +} // namespace PHARE::hdf5::h5 + + +#endif /* _PHARE_HDF5_DETAIL_H5_GROUP_SCANNER_HPP_ */ diff --git a/src/hdf5/detail/h5/h5_file.hpp b/src/hdf5/detail/h5/h5_file.hpp index 06ec1ec55..a72124fa4 100644 --- a/src/hdf5/detail/h5/h5_file.hpp +++ b/src/hdf5/detail/h5/h5_file.hpp @@ -10,6 +10,11 @@ #include "core/utilities/mpi_utils.hpp" #include "core/utilities/meta/meta_utilities.hpp" +#include + +#include "group_scanner.hpp" + + namespace PHARE::hdf5::h5 { using HiFile = HighFive::File; @@ -39,6 +44,7 @@ NO_DISCARD auto vector_for_dim() return std::vector>>(); } + class HighFiveFile { public: @@ -69,7 +75,8 @@ class HighFiveFile ~HighFiveFile() {} - NO_DISCARD HiFile& file() { return h5file_; } + NO_DISCARD auto& file() { return h5file_; } + NO_DISCARD auto& file() const { return h5file_; } template @@ -245,6 +252,13 @@ class HighFiveFile } + + std::unordered_set scan_for_groups(std::vector const& contains) + { + return GroupScanner{*this, contains}.scan(); + } + + HighFiveFile(const HighFiveFile&) = delete; HighFiveFile(const HighFiveFile&&) = delete; HighFiveFile& operator=(const HighFiveFile&) = delete; @@ -287,7 +301,6 @@ class HighFiveFile - } // namespace PHARE::hdf5::h5 diff --git a/src/phare_core.hpp b/src/phare_core.hpp index 0d6ef0d83..2e04f7686 100644 --- a/src/phare_core.hpp +++ b/src/phare_core.hpp @@ -55,9 +55,6 @@ struct PHARE_Types = PHARE::core::IonPopulation; using Ions_t = PHARE::core::Ions; using Electrons_t = PHARE::core::Electrons; - - using ParticleInitializerFactory - = PHARE::core::ParticleInitializerFactory; }; struct PHARE_Sim_Types diff --git a/src/simulator/phare_types.hpp b/src/simulator/phare_types.hpp index 08d2c7707..b336d1697 100644 --- a/src/simulator/phare_types.hpp +++ b/src/simulator/phare_types.hpp @@ -29,9 +29,9 @@ struct PHARE_Types using MaxwellianParticleInitializer_t = typename core_types::MaxwellianParticleInitializer_t; using IonPopulation_t = typename core_types::IonPopulation_t; using Electrons_t = typename core_types::Electrons_t; - using ParticleInitializerFactory = typename core_types::ParticleInitializerFactory; - + using ParticleInitializerFactory + = amr::ParticleInitializerFactory; using amr_types = PHARE::amr::PHARE_Types; using hierarchy_t = typename amr_types::hierarchy_t; @@ -40,7 +40,6 @@ struct PHARE_Types - using solver_types = PHARE::solver::PHARE_Types; using IPhysicalModel = typename solver_types::IPhysicalModel; using HybridModel_t = typename solver_types::HybridModel_t; diff --git a/tests/core/data/particle_initializer/CMakeLists.txt b/tests/amr/data/particles/initializer/CMakeLists.txt similarity index 96% rename from tests/core/data/particle_initializer/CMakeLists.txt rename to tests/amr/data/particles/initializer/CMakeLists.txt index a380b41b7..fc7506044 100644 --- a/tests/core/data/particle_initializer/CMakeLists.txt +++ b/tests/amr/data/particles/initializer/CMakeLists.txt @@ -11,7 +11,7 @@ target_include_directories(${PROJECT_NAME} PRIVATE ) target_link_libraries(${PROJECT_NAME} PRIVATE - phare_core + phare_amr phare_initializer ${GTEST_LIBS}) diff --git a/tests/core/data/particle_initializer/test_main.cpp b/tests/amr/data/particles/initializer/test_main.cpp similarity index 77% rename from tests/core/data/particle_initializer/test_main.cpp rename to tests/amr/data/particles/initializer/test_main.cpp index 405d42e5a..4008b57e1 100644 --- a/tests/core/data/particle_initializer/test_main.cpp +++ b/tests/amr/data/particles/initializer/test_main.cpp @@ -4,10 +4,10 @@ #include "core/utilities/span.hpp" #include "core/data/grid/gridlayout.hpp" #include "core/data/grid/gridlayoutimplyee.hpp" -#include "core/data/ions/particle_initializers/particle_initializer_factory.hpp" #include "core/data/particles/particle_array.hpp" #include "initializer/data_provider.hpp" +#include "amr/data/particles/initializers/particle_initializer_factory.hpp" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -18,9 +18,6 @@ using namespace PHARE::core; using namespace PHARE::initializer; -using namespace PHARE::core; -using namespace PHARE::initializer; - #include "tests/initializer/init_functions.hpp" using namespace PHARE::initializer::test_fn::func_1d; // density/etc are here @@ -32,7 +29,7 @@ using ParticleArrayT = ParticleArray<1>; TEST(AParticleIinitializerFactory, takesAPHAREDictToCreateAParticleVectorInitializer) { PHARE::initializer::PHAREDict dict; - dict["name"] = std::string{"MaxwellianParticleInitializer"}; + dict["name"] = std::string{"maxwellian"}; dict["fn_type"] = std::string{"vector"}; dict["density"] = static_cast>(density); dict["bulk_velocity_x"] = static_cast>(vx); @@ -43,11 +40,12 @@ TEST(AParticleIinitializerFactory, takesAPHAREDictToCreateAParticleVectorInitial dict["thermal_velocity_y"] = static_cast>(vthy); dict["thermal_velocity_z"] = static_cast>(vthz); - dict["charge"] = 1.; - dict["nbrPartPerCell"] = int{100}; - dict["basis"] = std::string{"Cartesian"}; + dict["charge"] = 1.; + dict["nbr_part_per_cell"] = int{100}; + dict["basis"] = std::string{"cartesian"}; - auto initializer = ParticleInitializerFactory::create(dict); + auto initializer + = PHARE::amr::ParticleInitializerFactory::create(dict); } int main(int argc, char** argv) diff --git a/tests/amr/models/test_models.cpp b/tests/amr/models/test_models.cpp index 7c902de4e..0dc83537b 100644 --- a/tests/amr/models/test_models.cpp +++ b/tests/amr/models/test_models.cpp @@ -58,11 +58,10 @@ using InitFunctionT = PHARE::initializer::InitFunction<1>; PHARE::initializer::PHAREDict createDict() { PHARE::initializer::PHAREDict dict; - dict["ions"]["nbrPopulations"] = std::size_t{2}; - dict["ions"]["pop0"]["name"] = std::string{"protons"}; - dict["ions"]["pop0"]["mass"] = 1.; - dict["ions"]["pop0"]["particle_initializer"]["name"] - = std::string{"MaxwellianParticleInitializer"}; + dict["ions"]["nbrPopulations"] = std::size_t{2}; + dict["ions"]["pop0"]["name"] = std::string{"protons"}; + dict["ions"]["pop0"]["mass"] = 1.; + dict["ions"]["pop0"]["particle_initializer"]["name"] = std::string{"maxwellian"}; dict["ions"]["pop0"]["particle_initializer"]["density"] = static_cast(density); dict["ions"]["pop0"]["particle_initializer"]["bulk_velocity_x"] @@ -85,14 +84,13 @@ PHARE::initializer::PHAREDict createDict() = static_cast(vthz); - dict["ions"]["pop0"]["particle_initializer"]["nbrPartPerCell"] = int{100}; - dict["ions"]["pop0"]["particle_initializer"]["charge"] = -1.; - dict["ions"]["pop0"]["particle_initializer"]["basis"] = std::string{"Cartesian"}; + dict["ions"]["pop0"]["particle_initializer"]["nbr_part_per_cell"] = int{100}; + dict["ions"]["pop0"]["particle_initializer"]["charge"] = -1.; + dict["ions"]["pop0"]["particle_initializer"]["basis"] = std::string{"cartesian"}; - dict["ions"]["pop1"]["name"] = std::string{"alpha"}; - dict["ions"]["pop1"]["mass"] = 1.; - dict["ions"]["pop1"]["particle_initializer"]["name"] - = std::string{"MaxwellianParticleInitializer"}; + dict["ions"]["pop1"]["name"] = std::string{"alpha"}; + dict["ions"]["pop1"]["mass"] = 1.; + dict["ions"]["pop1"]["particle_initializer"]["name"] = std::string{"maxwellian"}; dict["ions"]["pop1"]["particle_initializer"]["density"] = static_cast(density); dict["ions"]["pop1"]["particle_initializer"]["bulk_velocity_x"] @@ -115,9 +113,9 @@ PHARE::initializer::PHAREDict createDict() = 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"}; + dict["ions"]["pop1"]["particle_initializer"]["nbr_part_per_cell"] = int{100}; + dict["ions"]["pop1"]["particle_initializer"]["charge"] = -1.; + dict["ions"]["pop1"]["particle_initializer"]["basis"] = std::string{"cartesian"}; dict["electromag"]["name"] = std::string{"EM"}; dict["electromag"]["electric"]["name"] = std::string{"E"}; diff --git a/tests/amr/tagging/test_tagging.cpp b/tests/amr/tagging/test_tagging.cpp index 238126749..6096e0f8c 100644 --- a/tests/amr/tagging/test_tagging.cpp +++ b/tests/amr/tagging/test_tagging.cpp @@ -179,9 +179,18 @@ struct TestTagger : public ::testing::Test { using gridlayout_type = GridLayout>; static auto constexpr dimension = dim; - HybridState state; + + PHARE::initializer::PHAREDict dict_; + HybridState state{dict_}; + + void init(Electromag& em, GridLayoutT const& layout) + { + ElectromagUserFuncInitializer{dict_["electromag"]}.init( + em, layout); + } }; + PHARE::initializer::PHAREDict dict_; GridLayoutT layout; UsableVecField B, E; @@ -190,15 +199,17 @@ struct TestTagger : public ::testing::Test std::vector tags; TestTagger() - : layout{TestGridLayout::make(20)} + : dict_{createDict()} + , layout{TestGridLayout::make(20)} , B{"EM_B", layout, HybridQuantity::Vector::B} , E{"EM_E", layout, HybridQuantity::Vector::E} - , model{createDict()} + , model{dict_} , tags(20 + layout.nbrGhosts(PHARE::core::QtyCentering::dual)) { B.set_on(model.state.electromag.B); E.set_on(model.state.electromag.E); - model.state.electromag.initialize(layout); + + model.init(model.state.electromag, layout); } }; diff --git a/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp b/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp index f3f4899cb..b954b4622 100644 --- a/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp +++ b/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp @@ -51,7 +51,7 @@ TEST_F(AMaxwellianParticleInitializer1D, loadsTheCorrectNbrOfParticles) { auto nbrCells = layout.nbrCells(); auto expectedNbrParticles = nbrParticlesPerCell * nbrCells[0]; - initializer->loadParticles(particles, layout); + initializer->loadParticles(particles, layout, "protons"); EXPECT_EQ(expectedNbrParticles, particles.size()); } @@ -60,7 +60,7 @@ TEST_F(AMaxwellianParticleInitializer1D, loadsTheCorrectNbrOfParticles) TEST_F(AMaxwellianParticleInitializer1D, loadsParticlesInTheDomain) { - initializer->loadParticles(particles, layout); + initializer->loadParticles(particles, layout, "protons"); for (auto const& particle : particles) { EXPECT_TRUE(particle.iCell[0] >= 50 && particle.iCell[0] <= 99); diff --git a/tests/core/numerics/ion_updater/CMakeLists.txt b/tests/core/numerics/ion_updater/CMakeLists.txt index e206a3e58..bdbb17a97 100644 --- a/tests/core/numerics/ion_updater/CMakeLists.txt +++ b/tests/core/numerics/ion_updater/CMakeLists.txt @@ -11,7 +11,7 @@ target_include_directories(${PROJECT_NAME} PRIVATE ) target_link_libraries(${PROJECT_NAME} PRIVATE - phare_core + phare_amr phare_simulator ${GTEST_LIBS}) diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 2d939cf7d..29efe7f20 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -3,6 +3,8 @@ #include "phare_core.hpp" #include "core/numerics/ion_updater/ion_updater.hpp" +#include "amr/data/electromag/electromag_initializer.hpp" +#include "amr/data/particles/initializers/particle_initializer_factory.hpp" #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" #include "tests/core/data/tensorfield/test_tensorfield_fixtures.hpp" @@ -203,13 +205,14 @@ struct ElectromagBuffers template struct IonsBuffers { - using PHARETypes = PHARE::core::PHARE_Types; - using UsableVecFieldND = UsableVecField; - using Grid = typename PHARETypes::Grid_t; - using GridLayout = typename PHARETypes::GridLayout_t; - using Ions = typename PHARETypes::Ions_t; - using ParticleArray = typename PHARETypes::ParticleArray_t; - using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory; + using PHARETypes = PHARE::core::PHARE_Types; + using UsableVecFieldND = UsableVecField; + using Grid = typename PHARETypes::Grid_t; + using GridLayout = typename PHARETypes::GridLayout_t; + using Ions = typename PHARETypes::Ions_t; + using ParticleArray = typename PHARETypes::ParticleArray_t; + using ParticleInitializerFactory + = PHARE::amr::ParticleInitializerFactory; Grid ionDensity; Grid ionMassDensity; @@ -352,8 +355,10 @@ struct IonUpdaterTest : public ::testing::Test using Electromag = typename PHARETypes::Electromag_t; using GridLayout = typename PHARE::core::GridLayout>; using ParticleArray = typename PHARETypes::ParticleArray_t; - using ParticleInitializerFactory = typename PHARETypes::ParticleInitializerFactory; + using ParticleInitializerFactory + = PHARE::amr::ParticleInitializerFactory; + using ElectromagInitializerFactory_t = PHARE::amr::ElectromagInitializerFactory; using IonUpdater = typename PHARE::core::IonUpdater; @@ -391,13 +396,13 @@ struct IonUpdaterTest : public ::testing::Test // now let's initialize Electromag fields to user input functions // and ion population particles to user supplied moments - - EM.initialize(layout); + ElectromagInitializerFactory_t::create(init_dict["electromag"]) + ->init(EM, layout); for (auto& pop : ions) { auto const& info = pop.particleInitializerInfo(); auto particleInitializer = ParticleInitializerFactory::create(info); - particleInitializer->loadParticles(pop.domainParticles(), layout); + particleInitializer->loadParticles(pop.domainParticles(), layout, pop.name()); } @@ -515,7 +520,7 @@ struct IonUpdaterTest : public ::testing::Test } } // end 1D - } // end pop loop + } // end pop loop PHARE::core::depositParticles(ions, layout, Interpolator{}, PHARE::core::DomainDeposit{}); diff --git a/tests/simulator/CMakeLists.txt b/tests/simulator/CMakeLists.txt index 86274d16f..689fdf1ea 100644 --- a/tests/simulator/CMakeLists.txt +++ b/tests/simulator/CMakeLists.txt @@ -14,10 +14,10 @@ add_no_mpi_python3_test(periodicity test_init_periodicity.py ${CMAKE_CURRENT_BIN if(HighFive) ## These test use dump diagnostics so require HighFive! - phare_python3_exec(9 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 - phare_python3_exec(9 restarts test_restarts.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 - phare_python3_exec(9 run test_run.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 - # phare_python3_exec(9 tagging test_tagging.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 + phare_python3_exec(9 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 + phare_python3_exec(9 restarts test_restarts.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 + phare_python3_exec(9 run test_run.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 + # phare_python3_exec(9 tagging test_tagging.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 if(testMPI) phare_mpi_python3_exec(9 3 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) @@ -25,9 +25,13 @@ if(HighFive) # doesn't make sense in serial phare_mpi_python3_exec(9 3 load_balancing test_load_balancing.py ${CMAKE_CURRENT_BINARY_DIR}) + phare_mpi_python3_exec(9 3 initing test_init_from_restart.py ${CMAKE_CURRENT_BINARY_DIR}) + endif(testMPI) phare_python3_exec(11, test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR}) + + endif() configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config.py ${CMAKE_CURRENT_BINARY_DIR}/config.py @ONLY) @@ -35,4 +39,3 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/config.py ${CMAKE_CURRENT_BINARY_DIR} add_subdirectory(initialize) add_subdirectory(refinement) add_subdirectory(advance) - diff --git a/tests/simulator/test_init_from_restart.py b/tests/simulator/test_init_from_restart.py new file mode 100644 index 000000000..68b71fff7 --- /dev/null +++ b/tests/simulator/test_init_from_restart.py @@ -0,0 +1,84 @@ +import sys +import copy +import unittest +import subprocess +import numpy as np +import pyphare.pharein as ph +from pathlib import Path + +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.hierarchy.patchdata import FieldData, ParticleData +from pyphare.pharesee.hierarchy.fromh5 import get_all_available_quantities_from_h5 +from pyphare.pharesee.hierarchy.hierarchy import format_timestamp +from pyphare.pharesee.hierarchy.hierarchy_utils import single_patch_for_LO +from pyphare.pharesee.hierarchy.hierarchy_utils import hierarchy_compare +from tests.simulator import SimulatorTest, test_restarts +from tests.diagnostic import dump_all_diags + + +time_step = 0.001 +time_step_nbr = 5 +final_time = time_step_nbr * time_step +first_mpi_size = 5 +ppc = 100 +cells = 200 +first_out = "test_init_from_restart" +secnd_out = "phare_outputs/reinit/secnd" +timestamps = np.arange(0, final_time + time_step, time_step) +restart_idx = Z = 2 +simInitArgs = dict( + time_step_nbr=time_step_nbr, + time_step=time_step, + cells=cells, + dl=0.3, + init_options=dict(dir=f"{first_out}/00000.00{Z}00", mpi_size=first_mpi_size), + diag_options=dict(format="phareh5", options=dict(dir=secnd_out, mode="overwrite")), +) + + +def setup_model(sim): + model = ph.MaxwellianFluidModel( + protons={"mass": 1, "charge": 1, "nbr_part_per_cell": ppc}, + alpha={"mass": 4, "charge": 1, "nbr_part_per_cell": ppc}, + ) + ph.ElectronModel(closure="isothermal", Te=0.12) + dump_all_diags(model.populations, timestamps=timestamps) + return model + + +class RestartsParserTest(SimulatorTest): + def __init__(self, *args, **kwargs): + super(RestartsParserTest, self).__init__(*args, **kwargs) + self.simulator = None + + def tearDown(self): + super(RestartsParserTest, self).tearDown() + if self.simulator is not None: + self.simulator.reset() + self.simulator = None + ph.global_vars.sim = None + + def test_reinit(self): + self.register_diag_dir_for_cleanup("phare_outputs/reinit") + sim = ph.Simulation(**copy.deepcopy(simInitArgs)) + setup_model(sim) + Simulator(sim).run().reset() + fidx, sidx = 4, 2 + datahier0 = get_all_available_quantities_from_h5(first_out, timestamps[fidx]) + datahier0.time_hier = { # swap times + format_timestamp(timestamps[sidx]): datahier0.time_hier[ + format_timestamp(timestamps[fidx]) + ] + } + datahier1 = get_all_available_quantities_from_h5(secnd_out, timestamps[sidx]) + qties = None + skip = ["protons_patchGhost", "alpha_patchGhost"] + ds = [single_patch_for_LO(d, qties, skip) for d in [datahier0, datahier1]] + self.assertTrue(hierarchy_compare(*ds, atol=1e-12)) + + +if __name__ == "__main__": + if Path("test_init_from_restart").exists(): + unittest.main() + else: + print("No source data found - see tools/test_data_gen.py ") diff --git a/tools/test_data_gen.py b/tools/test_data_gen.py new file mode 100644 index 000000000..3be088070 --- /dev/null +++ b/tools/test_data_gen.py @@ -0,0 +1,48 @@ +import sys +import copy +import unittest +import subprocess +import numpy as np +import pyphare.pharein as ph +from pyphare.simulator.simulator import Simulator +from tests.simulator import test_restarts +from tests.diagnostic import dump_all_diags + +output_dir = "phare_data_gen/" + + +def run_restartable_data_sim(): + from tests.simulator import test_init_from_restart as restartable_sim + + """uses params from tests_restarts.py""" + out = output_dir + "/tests/simulator/test_init_from_restart" + simput = copy.deepcopy(test_restarts.simArgs) + simput["restart_options"]["dir"] = out + simput["restart_options"]["timestamps"] = restartable_sim.timestamps + simput["diag_options"]["options"]["dir"] = out + sim = ph.Simulation(**simput) + dump_all_diags( + test_restarts.setup_model().populations, timestamps=restartable_sim.timestamps + ) + Simulator(sim).run() + + +def launch(fn, n=5): + """Launch secondary process to run first simulation to avoid initalizing MPI now""" + cmd = f"mpirun -n {n} python3 -O {__file__} {fn}" + print(cmd) + try: + p = subprocess.run(cmd.split(" "), check=True, capture_output=True) + except subprocess.CalledProcessError as e: + print("CalledProcessError", e) + + +if __name__ == "__main__": + if len(sys.argv) > 1: + globals()[sys.argv[1]]() + else: + for k in [k for k in list(globals().keys()) if k.startswith("run_")]: + launch(k) +# +# tar czf phare_data_gen.tar.xz phare_data_gen +# mc put phare_data_gen.tar.xz minio/phare/phare_data_gen.tar.xz --insecure