wip

Author: rochSmets

.gdbinit

@@ -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")

pyphare/pyphare/pharesee/hierarchy/__init__.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}")

pyphare/pyphare/pharesee/hierarchy/fromfunc.py

@@ -36,7 +36,7 @@ void DefaultHybridTaggerStrategy<HybridModel>::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, _]

src/amr/tagging/default_hybrid_tagger_strategy.hpp

@@ -925,13 +925,15 @@ 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);
             auto& fy = flux.getComponent(Component::Y);
             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)));

tests/core/numerics/ion_updater/test_updater.cpp

@@ -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()