diff --git a/.gdbinit b/.gdbinit deleted file mode 100644 index 9c8c9d7a1..000000000 --- a/.gdbinit +++ /dev/null @@ -1,3 +0,0 @@ -break __assert_fail -break abort -catch throw diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py index 318b1a24c..280428513 100644 --- a/pyphare/pyphare/pharesee/hierarchy/__init__.py +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -41,6 +41,8 @@ def hierarchy_from( return hierarchy_from_sim(simulator, qty, pop=pop) if func is not None and hier is not None: + if not callable(func): + raise TypeError("func must be callable") return hierarchy_from_func(func, hier, **kwargs) raise ValueError("can't make hierarchy") diff --git a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py index 5d270d37f..ac313c298 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py @@ -1,7 +1,5 @@ from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from -import numpy as np - def hierarchy_from_func1d(func, hier, **kwargs): assert hier.ndim == 1 @@ -33,8 +31,33 @@ def compute_(patch_datas, **kwargs): def hierarchy_from_func(func, hier, **kwargs): + + """ + Route hierarchical computation to appropriate dimension handler. + + Parameters + ---------- + func : callable + Function to apply to coordinates of the hierarchy + hier : Hierarchy + Hierarchy object (1D or 2D) + **kwargs : dict + Additional arguments passed to func + + Returns + ------- + dict + Computed hierarchical data + + Raises + ------ + ValueError + If hierarchy dimension is not supported + """ + if hier.ndim == 1: return hierarchy_from_func1d(func, hier, **kwargs) if hier.ndim == 2: return hierarchy_from_func2d(func, hier, **kwargs) - + else: + raise ValueError(f"Unsupported hierarchy dimension: {hier.ndim}") diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp index 741a52b4e..1bf22f723 100644 --- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp +++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp @@ -36,7 +36,7 @@ void DefaultHybridTaggerStrategy::tag(HybridModel& model, auto& By = model.state.electromag.B.getComponent(PHARE::core::Component::Y); auto& Bz = model.state.electromag.B.getComponent(PHARE::core::Component::Z); - auto& N = model.state.ions.chargeDensity(); + // auto& N = model.state.ions.chargeDensity(); // we loop on cell indexes for all qties regardless of their centering auto const& [start_x, _] diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 09013e517..2ab71bfb4 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -925,6 +925,7 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments) for (auto ix = ix0; ix <= ix1; ++ix) { auto& density = pop.particleDensity(); + auto& chargeDensity = pop.chargeDensity(); auto& flux = pop.flux(); auto& fx = flux.getComponent(Component::X); @@ -932,6 +933,7 @@ TYPED_TEST(IonUpdaterTest, thatNoNaNsExistOnPhysicalNodesMoments) auto& fz = flux.getComponent(Component::Z); EXPECT_FALSE(std::isnan(density(ix))); + EXPECT_FALSE(std::isnan(chargeDensity(ix))); EXPECT_FALSE(std::isnan(fx(ix))); EXPECT_FALSE(std::isnan(fy(ix))); EXPECT_FALSE(std::isnan(fz(ix))); diff --git a/tests/simulator/initialize/test_density_init.py b/tests/simulator/initialize/test_density_init.py index b7de9850f..a284c8756 100644 --- a/tests/simulator/initialize/test_density_init.py +++ b/tests/simulator/initialize/test_density_init.py @@ -24,8 +24,9 @@ def ions_mass_density_func1d(x, **kwargs): masses = kwargs["masses"] # list of float : the ion pop masses densities = kwargs["densities"] # list of callable : the ion pop density profiles + if len(masses) != len(densities): + raise ValueError("Length of masses and densities must be equal") - assert len(masses) == len(densities) funcs = np.zeros((x.size, len(masses))) for i, (mass, density) in enumerate(zip(masses, densities)): @@ -37,8 +38,8 @@ def ions_mass_density_func1d(x, **kwargs): def ions_charge_density_func1d(x, **kwargs): charges = kwargs["charges"] # list of float : the ion pop charges densities = kwargs["densities"] # list of callable : the ion pop density profiles - - assert len(charges) == len(densities) + if len(charges) != len(densities): + raise ValueError("Length of charges and densities must be equal") funcs = np.zeros((x.size, len(charges))) @@ -52,10 +53,11 @@ def ions_charge_density_func1d(x, **kwargs): def ions_mass_density_func2d(x, y, **kwargs): masses = kwargs["masses"] # list of float : the ion pop masses densities = kwargs["densities"] # list of callable : the ion pop density profiles + if len(masses) != len(densities): + raise ValueError("Length of masses and densities must be equal") yv, xv = np.meshgrid(y, x) - assert len(masses) == len(densities) funcs = np.zeros((x.size, y.size, len(masses))) for i, (mass, density) in enumerate(zip(masses, densities)): @@ -67,10 +69,11 @@ def ions_mass_density_func2d(x, y, **kwargs): def ions_charge_density_func2d(x, y, **kwargs): charges = kwargs["charges"] # list of float : the ion pop charges densities = kwargs["densities"] # list of callable : the ion pop density profiles + if len(charges) != len(densities): + raise ValueError("Length of charges and densities must be equal") yv, xv = np.meshgrid(y, x) - assert len(charges) == len(densities) funcs = np.zeros((x.size, y.size, len(charges))) for i, (charge, density) in enumerate(zip(charges, densities)): @@ -140,8 +143,8 @@ def config_1d(): bx=bx_1d, by=by_1d, bz=bz_1d, - main={"mass": masses[0], "charge": charges[0], "density": densityMain_1d, "nbr_part_per_cell": 1000, **v_pop}, - beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_1d, "nbr_part_per_cell": 1000, **v_pop}, + main={"mass": masses[0], "charge": charges[0], "density": densityMain_1d, "nbr_part_per_cell": 100, **v_pop}, + beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_1d, "nbr_part_per_cell": 100, **v_pop}, ) ph.ElectronModel(closure="isothermal", Te=0.0) @@ -222,8 +225,8 @@ def config_2d(): bx=bx_2d, by=by_2d, bz=bz_2d, - main={"mass": masses[0], "charge": charges[0], "density": densityMain_2d, "nbr_part_per_cell": 2000, **v_pop}, - beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_2d, "nbr_part_per_cell": 2000, **v_pop}, + main={"mass": masses[0], "charge": charges[0], "density": densityMain_2d, "nbr_part_per_cell": 100, **v_pop}, + beam={"mass": masses[1], "charge": charges[1], "density": densityBeam_2d, "nbr_part_per_cell": 100, **v_pop}, ) ph.ElectronModel(closure="isothermal", Te=0.0) @@ -255,23 +258,21 @@ def main(): ph.global_vars.sim = None Simulator(config_2d()).run().reset() - - def assert_close_enough(h, H): + def noise_level(h, H): + noises = list() for lvl_h, lvl_H in zip(h.levels(time).values(), H.levels(time).values()): for patch_h, patch_H in zip(lvl_h.patches, lvl_H.patches): pd_h = patch_h.patch_datas["value"] pd_H = patch_H.patch_datas["value"] - ghosts_num = pd_h.ghosts_nbr[0] - if pd_H.ndim == 1: - dset_h = pd_h.dataset[ghosts_num:-ghosts_num] - dset_H = pd_H.dataset[ghosts_num:-ghosts_num] - if pd_H.ndim == 2: - dset_h = pd_h.dataset[ghosts_num:-ghosts_num, ghosts_num:-ghosts_num] - dset_H = pd_H.dataset[ghosts_num:-ghosts_num, ghosts_num:-ghosts_num] + dset_h = pd_h[pd_h.box] + dset_H = pd_H[pd_H.box] + + noise = np.std(dset_h - dset_H) + # print(f"noise = ", noise) + noises.append(noise) + return noises - for h_, H_ in zip(dset_h, dset_H): - np.testing.assert_almost_equal(h_, H_, decimal=1) fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(6, 8)) @@ -287,8 +288,8 @@ def assert_close_enough(h, H): H1 = hierarchy_from(hier=h1, func=ions_mass_density_func1d, masses=masses, densities=(densityMain_1d, densityBeam_1d)) H2 = hierarchy_from(hier=h2, func=ions_charge_density_func1d, charges=charges, densities=(densityMain_1d, densityBeam_1d)) - assert_close_enough(h1, H1) - assert_close_enough(h2, H2) + assert np.mean(noise_level(h1, H1)) < 0.18 + assert np.mean(noise_level(h2, H2)) < 0.12 cycle = plt.rcParams['axes.prop_cycle'].by_key()['color'] @@ -313,8 +314,8 @@ def assert_close_enough(h, H): H1 = hierarchy_from(hier=h1, func=ions_mass_density_func2d, masses=masses, densities=(densityMain_2d, densityBeam_2d)) H2 = hierarchy_from(hier=h2, func=ions_charge_density_func2d, charges=charges, densities=(densityMain_2d, densityBeam_2d)) - assert_close_enough(h1, H1) - assert_close_enough(h2, H2) + assert np.mean(noise_level(h1, H1)) < 0.18 + assert np.mean(noise_level(h2, H2)) < 0.12 cmap = mpl.colormaps['viridis'] @@ -327,8 +328,5 @@ def assert_close_enough(h, H): plt.savefig("nCheck.pdf", dpi=300) - - # /home/smets/codes/far/PHARE/tests/simulator/initialize - if __name__ == "__main__": main()