PR comments

pyphare/pyphare/pharein/diagnostics.py

@@ -222,6 +222,7 @@ def population_in_model(population):
 class FluidDiagnostics_(Diagnostics):
     fluid_quantities = [
         "density",
+        "charge_density",
         "mass_density",
         "flux",
         "bulkVelocity",

pyphare/pyphare/pharesee/hierarchy/__init__.py

@@ -11,10 +11,11 @@
 
 
 def hierarchy_from(
-    simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs
+    simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, func=None, **kwargs
 ):
     from .fromh5 import hierarchy_fromh5
     from .fromsim import hierarchy_from_sim
+    from .fromfunc import hierarchy_from_func
 
     """
     this function reads an HDF5 PHARE file and returns a PatchHierarchy from
@@ -39,4 +40,9 @@ def hierarchy_from(
     if simulator is not None and qty is not None:
         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/fromfunc.py

@@ -0,0 +1,63 @@
+from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+
+def hierarchy_from_func1d(func, hier, **kwargs):
+    assert hier.ndim == 1
+
+    def compute_(patch_datas, **kwargs):
+        ref_name = next(iter(patch_datas.keys()))
+        x_ = patch_datas[ref_name].x
+
+        return (
+            {"name": "value", "data": func(x_, **kwargs), "centering": patch_datas[ref_name].centerings},
+        )
+
+    return compute_hier_from(compute_, hier, **kwargs)
+
+
+def hierarchy_from_func2d(func, hier, **kwargs):
+    from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+    def compute_(patch_datas, **kwargs):
+        ref_name = next(iter(patch_datas.keys()))
+        x_ = patch_datas[ref_name].x
+        y_ = patch_datas[ref_name].y
+
+        return (
+            {"name": "value", "data": func(x_, y_, **kwargs), "centering": patch_datas[ref_name].centerings},
+        )
+
+    return compute_hier_from(compute_, hier, **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/hierarchy.py

@@ -430,8 +430,7 @@ def plot1d(self, **kwargs):
                 if qty is None:
                     qty = pdata_names[0]
 
-                layout = patch.patch_datas[qty].layout
-                nbrGhosts = layout.nbrGhostFor(qty)
+                nbrGhosts = patch.patch_datas[qty].ghosts_nbr
                 val = patch.patch_datas[qty][patch.box]
                 x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]]
                 label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip)
@@ -542,9 +541,10 @@ def plot2d(self, **kwargs):
         if "ylim" in kwargs:
             ax.set_ylim(kwargs["ylim"])
 
-        divider = make_axes_locatable(ax)
-        cax = divider.append_axes("right", size="5%", pad=0.08)
-        plt.colorbar(im, ax=ax, cax=cax)
+        if kwargs.get("cbar", True):
+            divider = make_axes_locatable(ax)
+            cax = divider.append_axes("right", size="5%", pad=0.08)
+            fig.colorbar(im, ax=ax, cax=cax)
 
         if kwargs.get("legend", None) is not None:
             ax.legend()

pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

@@ -36,6 +36,7 @@
     "momentum_tensor_zz": "Mzz",
     "density": "rho",
     "mass_density": "rho",
+    "charge_density": "rho",
     "tags": "tags",
 }
 

pyphare/pyphare/pharesee/plotting.py

@@ -274,7 +274,7 @@ def finest_field_plot(run_path, qty, **kwargs):
             time = times[0]
         interpolator, finest_coords = r.GetVi(time, merged=True, interp=interp)[qty]
     elif qty == "rho":
-        file = os.path.join(run_path, "ions_density.h5")
+        file = os.path.join(run_path, "ions_charge_density.h5")
         if time is None:
             times = get_times_from_h5(file)
             time = times[0]

pyphare/pyphare/pharesee/run/run.py

@@ -27,7 +27,7 @@
     "EM_B": "B",
     "EM_E": "E",
     "ions_bulkVelocity": "Vi",
-    "ions_density": "Ni",
+    "ions_charge_density": "Ni",
     "particle_count": "nppc",
 }
 
@@ -114,7 +114,7 @@ def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs):
         return ScalarField(self._get(hier, time, merged, interp))
 
     def GetNi(self, time, merged=False, interp="nearest", **kwargs):
-        hier = self._get_hierarchy(time, "ions_density.h5", **kwargs)
+        hier = self._get_hierarchy(time, "ions_charge_density.h5", **kwargs)
         return ScalarField(self._get(hier, time, merged, interp))
 
     def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs):
@@ -151,7 +151,7 @@ def GetPi(self, time, merged=False, interp="nearest", **kwargs):
         return self._get(Pi, time, merged, interp)  # should later be a TensorField
 
     def GetPe(self, time, merged=False, interp="nearest", all_primal=True):
-        hier = self._get_hierarchy(time, "ions_density.h5")
+        hier = self._get_hierarchy(time, "ions_charge_density.h5")
 
         Te = hier.sim.electrons.closure.Te
 

pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py

@@ -177,11 +177,11 @@ def test_large_patchoverlaps(self, expected):
             collections=[
                 {
                     "boxes": overlap_boxes,
-                    "facecolor": "yellow",
+                    "facecolor": "grey",
                 },
                 {
                     "boxes": [p.box for p in hierarchy.level(ilvl).patches],
-                    "facecolor": "grey",
+                    "facecolor": "yellow",
                 },
             ],
         )

pyphare/pyphare_tests/test_pharesee/test_hierarchy.py

@@ -145,7 +145,7 @@ def vthz(x, y):
             for quantity in ["E", "B"]:
                 ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
-            for quantity in ["density", "bulkVelocity"]:
+            for quantity in ["charge_density", "bulkVelocity"]:
                 ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps)
 
             return sim

src/amr/level_initializer/hybrid_level_initializer.hpp

@@ -98,7 +98,7 @@ namespace solver
                     core::depositParticles(ions, layout, interpolate_, core::LevelGhostDeposit{});
                 }
 
-                ions.computeDensity();
+                ions.computeChargeDensity();
                 ions.computeBulkVelocity();
             }
 

src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp

@@ -385,20 +385,21 @@ namespace amr
                 {
                     // first thing to do is to project patchGhostParitcles moments
                     auto& patchGhosts = pop.patchGhostParticles();
-                    auto& density     = pop.density();
+                    auto& particleDensity     = pop.particleDensity();
+                    auto& chargeDensity     = pop.chargeDensity();
                     auto& flux        = pop.flux();
 
-                    interpolate_(makeRange(patchGhosts), density, flux, layout);
+                    interpolate_(makeRange(patchGhosts), particleDensity, chargeDensity, flux, layout);
 
                     if (level.getLevelNumber() > 0) // no levelGhost on root level
                     {
                         // then grab levelGhostParticlesOld and levelGhostParticlesNew
                         // and project them with alpha and (1-alpha) coefs, respectively
                         auto& levelGhostOld = pop.levelGhostParticlesOld();
-                        interpolate_(makeRange(levelGhostOld), density, flux, layout, 1. - alpha);
+                        interpolate_(makeRange(levelGhostOld), particleDensity, chargeDensity, flux, layout, 1. - alpha);
 
                         auto& levelGhostNew = pop.levelGhostParticlesNew();
-                        interpolate_(makeRange(levelGhostNew), density, flux, layout, alpha);
+                        interpolate_(makeRange(levelGhostNew), particleDensity, chargeDensity, flux, layout, alpha);
                     }
                 }
             }
@@ -538,7 +539,7 @@ namespace amr
 
                 auto& J  = hybridModel.state.J;
                 auto& Vi = hybridModel.state.ions.velocity();
-                auto& Ni = hybridModel.state.ions.density();
+                auto& Ni = hybridModel.state.ions.chargeDensity();
 
                 Jold_.copyData(J);
                 ViOld_.copyData(Vi);

src/amr/physical_models/hybrid_model.hpp

@@ -154,7 +154,7 @@ void HybridModel<GridLayoutT, Electromag, Ions, Electrons, AMR_Types, Grid_t>::f
 
     hybridInfo.modelMagnetic        = core::VecFieldNames{state.electromag.B};
     hybridInfo.modelElectric        = core::VecFieldNames{state.electromag.E};
-    hybridInfo.modelIonDensity      = state.ions.densityName();
+    hybridInfo.modelIonDensity      = state.ions.chargeDensityName();
     hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()};
     hybridInfo.modelCurrent         = core::VecFieldNames{state.J};
 

src/amr/tagging/default_hybrid_tagger_strategy.hpp

@@ -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.density();
+    // auto& N = model.state.ions.chargeDensity();
 
     // we loop on cell indexes for all qties regardless of their centering
     auto const& [start_x, _]

src/core/data/electrons/electrons.hpp

@@ -65,7 +65,7 @@ class StandardHybridElectronFluxComputer
     {
         if (isUsable())
         {
-            return ions_.density();
+            return ions_.chargeDensity();
         }
         else
         {
@@ -78,7 +78,7 @@ class StandardHybridElectronFluxComputer
     {
         if (isUsable())
         {
-            return ions_.density();
+            return ions_.chargeDensity();
         }
         else
         {
@@ -112,23 +112,23 @@ class StandardHybridElectronFluxComputer
         auto const& Vix = ions_.velocity()(Component::X);
         auto const& Viy = ions_.velocity()(Component::Y);
         auto const& Viz = ions_.velocity()(Component::Z);
-        auto const& Ni  = ions_.density();
+        auto const& Ne  = ions_.chargeDensity();
 
         auto& Vex = Ve_(Component::X);
         auto& Vey = Ve_(Component::Y);
         auto& Vez = Ve_(Component::Z);
 
         // from Ni because all components defined on primal
-        layout.evalOnBox(Ni, [&](auto const&... args) {
+        layout.evalOnBox(Ne, [&](auto const&... args) {
             auto arr = std::array{args...};
 
             auto const JxOnVx = GridLayout::project(Jx, arr, GridLayout::JxToMoments());
             auto const JyOnVy = GridLayout::project(Jy, arr, GridLayout::JyToMoments());
             auto const JzOnVz = GridLayout::project(Jz, arr, GridLayout::JzToMoments());
 
-            Vex(arr) = Vix(arr) - JxOnVx / Ni(arr);
-            Vey(arr) = Viy(arr) - JyOnVy / Ni(arr);
-            Vez(arr) = Viz(arr) - JzOnVz / Ni(arr);
+            Vex(arr) = Vix(arr) - JxOnVx / Ne(arr);
+            Vey(arr) = Viy(arr) - JyOnVy / Ne(arr);
+            Vez(arr) = Viz(arr) - JzOnVz / Ne(arr);
         });
     }
 
@@ -211,7 +211,7 @@ class IsothermalElectronPressureClosure
         if (!Pe_.isUsable())
             throw std::runtime_error("Error - isothermal closure pressure not usable");
 
-        auto const& Ne_ = ions_.density();
+        auto const& Ne_ = ions_.chargeDensity();
         std::transform(std::begin(Ne_), std::end(Ne_), std::begin(Pe_),
                        [this](auto n) { return n * Te_; });
     }

src/core/data/ions/ion_population/ion_population.hpp

@@ -31,7 +31,8 @@ namespace core
             , mass_{initializer["mass"].template to<double>()}
             , flux_{name_ + "_flux", HybridQuantity::Vector::V}
             , momentumTensor_{name_ + "_momentumTensor", HybridQuantity::Tensor::M}
-            , rho_{name_ + "_rho", HybridQuantity::Scalar::rho}
+            , particleDensity_{name_ + "_particleDensity", HybridQuantity::Scalar::rho}
+            , chargeDensity_{name_ + "_chargeDensity", HybridQuantity::Scalar::rho}
             , particles_{name_}
             , particleInitializerInfo_{initializer["particle_initializer"]}
         {
@@ -42,20 +43,16 @@ namespace core
 
         NO_DISCARD std::string const& name() const { return name_; }
 
-
         NO_DISCARD auto const& particleInitializerInfo() const { return particleInitializerInfo_; }
 
-
-
         NO_DISCARD bool isUsable() const
         {
-            return core::isUsable(particles_, rho_, flux_, momentumTensor_);
+            return core::isUsable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_);
         }
 
-
         NO_DISCARD bool isSettable() const
         {
-            return core::isSettable(particles_, rho_, flux_, momentumTensor_);
+            return core::isSettable(particles_, particleDensity_, chargeDensity_, flux_, momentumTensor_);
         }
 
         NO_DISCARD auto& domainParticles() const { return particles_.domainParticles(); }
@@ -79,9 +76,11 @@ namespace core
             return particles_.levelGhostParticlesNew();
         }
 
+        NO_DISCARD field_type const& particleDensity() const { return particleDensity_; }
+        NO_DISCARD field_type& particleDensity() { return particleDensity_; }
 
-        NO_DISCARD field_type const& density() const { return rho_; }
-        NO_DISCARD field_type& density() { return rho_; }
+        NO_DISCARD field_type const& chargeDensity() const { return chargeDensity_; }
+        NO_DISCARD field_type& chargeDensity() { return chargeDensity_; }
 
         NO_DISCARD VecField const& flux() const { return flux_; }
         NO_DISCARD VecField& flux() { return flux_; }
@@ -100,7 +99,8 @@ namespace core
 
         NO_DISCARD auto getCompileTimeResourcesViewList()
         {
-            return std::forward_as_tuple(flux_, momentumTensor_, rho_, particles_);
+            return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_,
+                                         particles_);
         }
 
 
@@ -124,7 +124,8 @@ namespace core
         double mass_;
         VecField flux_;
         TensorField momentumTensor_;
-        field_type rho_;
+        field_type particleDensity_;
+        field_type chargeDensity_;
         ParticlesPack<ParticleArray> particles_;
         initializer::PHAREDict const& particleInitializerInfo_;
     };