fix

pyphare/pyphare/pharein/diagnostics.py

@@ -291,7 +291,7 @@ def __init__(self, **kwargs):
             if for_total_ions(**kwargs):
                 needed_quantities = ["mass_density", "bulkVelocity", "momentum_tensor"]
             else:
-                needed_quantities = ["charge_density", "flux", "momentum_tensor"]
+                needed_quantities = ["density", "flux", "momentum_tensor"]
 
             for quantity in needed_quantities:
                 kwargs["quantity"] = quantity

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -353,6 +353,8 @@ def box_to_Rectangle(self, box):
         return Rectangle(box.lower, *box.shape)
 
     def plot_2d_patches(self, ilvl, collections, **kwargs):
+        from matplotlib.patches import Rectangle
+
         if isinstance(collections, list) and all(
             [isinstance(el, Box) for el in collections]
         ):
@@ -363,35 +365,48 @@ def plot_2d_patches(self, ilvl, collections, **kwargs):
         level_domain_box = self.level_domain_box(ilvl)
         mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max()
 
-        fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6)))
+        fig, ax = kwargs.get("subplot", plt.subplots(figsize=(16, 16)))
+
+        color = 1
+        i0, j0 = level_domain_box.lower
+        i1, j1 = level_domain_box.upper
+        ij = np.zeros((i1 - i0 + 1, j1 - j0 + 1)) + np.nan
+        ix = np.arange(i0, i1 + 1)
+        iy = np.arange(j0, j1 + 1)
 
         for collection in collections:
-            facecolor = collection.get("facecolor", "none")
-            edgecolor = collection.get("edgecolor", "purple")
-            alpha = collection.get("alpha", 1)
-            rects = [self.box_to_Rectangle(box) for box in collection["boxes"]]
-
-            ax.add_collection(
-                PatchCollection(
-                    rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor
-                )
-            )
+            facecolor = collection.get("facecolor", np.nan)
+            for box in collection["boxes"]:
+                if isinstance(box, Box):
+                    i0, j0 = box.lower
+                    i1, j1 = box.upper
+                    ij[i0 : i1 + 1, j0 : j1 + 1] = facecolor
+                else:
+                    ij[box] = facecolor
+
+        ax.pcolormesh(ix, iy, ij.T, edgecolors="k", cmap="jet")
+        ax.set_xticks(ix)
+        ax.set_yticks(iy)
+
+        for patch in self.level(ilvl).patches:
+            box = patch.box
+            r = Rectangle(box.lower - 0.5, *(box.upper + 0.5))
+
+            r.set_edgecolor("r")
+            r.set_facecolor("none")
+            r.set_linewidth(2)
+            ax.add_patch(r)
 
         if "title" in kwargs:
             from textwrap import wrap
 
             xfigsize = int(fig.get_size_inches()[0] * 10)  # 10 characters per inch
             ax.set_title("\n".join(wrap(kwargs["title"], xfigsize)))
 
-        major_ticks = np.arange(mi - 5, ma + 5 + 5, 5)
-        ax.set_xticks(major_ticks)
-        ax.set_yticks(major_ticks)
-
-        minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1)
-        ax.set_xticks(minor_ticks, minor=True)
-        ax.set_yticks(minor_ticks, minor=True)
-
-        ax.grid(which="both")
+        if "xlim" in kwargs:
+            ax.set_xlim(kwargs["xlim"])
+        if "ylim" in kwargs:
+            ax.set_ylim(kwargs["ylim"])
 
         return fig
 

pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

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

src/amr/level_initializer/hybrid_level_initializer.hpp

@@ -98,7 +98,7 @@ namespace solver
                     core::depositParticles(ions, layout, interpolate_, core::LevelGhostDeposit{});
                 }
 
-                ions.computeParticleDensity(); // TODO ouam : do we need here the charge density
+                ions.computeChargeDensity();
                 ions.computeBulkVelocity();
             }
 

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.particleDensityName();
+    hybridInfo.modelIonDensity      = state.ions.chargeDensityName();
     hybridInfo.modelIonBulkVelocity = core::VecFieldNames{state.ions.velocity()};
     hybridInfo.modelCurrent         = core::VecFieldNames{state.J};
 

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
         {
@@ -102,7 +102,7 @@ class StandardHybridElectronFluxComputer
     }
 
 
-    void computeParticleDensity() {} // TODO ouam : which is the same as the charge particle... except the sign
+    void computeDensity() {}
 
     void computeBulkVelocity(GridLayout const& layout)
     {
@@ -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_; });
     }
@@ -281,7 +281,7 @@ class ElectronMomentModel
 
 
 
-    void computeParticleDensity() { fluxComput_.computeParticleDensity(); }
+    void computeDensity() { fluxComput_.computeDensity(); }
     void computeBulkVelocity(GridLayout const& layout) { fluxComput_.computeBulkVelocity(layout); }
     void computePressure(GridLayout const& layout) { pressureClosure_.computePressure(layout); }
 
@@ -310,7 +310,7 @@ class Electrons : public LayoutHolder<typename Ions::gridlayout_type>
     {
         if (isUsable())
         {
-            momentModel_.computeParticleDensity();
+            momentModel_.computeDensity();
             momentModel_.computeBulkVelocity(layout);
             momentModel_.computePressure(layout);
         }

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

@@ -50,13 +50,14 @@ namespace core
 
         NO_DISCARD bool isUsable() const
         {
-            return particles_.isUsable() && particleDensity_.isUsable() && chargeDensity_.isUsable() && flux_.isUsable()
-                   && momentumTensor_.isUsable();
+            return particles_.isUsable() && particleDensity_.isUsable() && chargeDensity_.isUsable()
+                   && flux_.isUsable() && momentumTensor_.isUsable();
         }
 
         NO_DISCARD bool isSettable() const
         {
-            return particles_.isSettable() && particleDensity_.isSettable() && chargeDensity_.isSettable() && flux_.isSettable()
+            return particles_.isSettable() && particleDensity_.isSettable()
+                   && chargeDensity_.isSettable() && flux_.isSettable()
                    && momentumTensor_.isSettable();
         }
 
@@ -81,9 +82,6 @@ namespace core
             return particles_.levelGhostParticlesNew();
         }
 
-        NO_DISCARD field_type const& density() const { return particleDensity_; }  // TODO ouam : to remove
-        NO_DISCARD field_type& density() { return particleDensity_; }  // TODO ouam : to remove
-
         NO_DISCARD field_type const& particleDensity() const { return particleDensity_; }
         NO_DISCARD field_type& particleDensity() { return particleDensity_; }
 
@@ -107,7 +105,8 @@ namespace core
 
         NO_DISCARD auto getCompileTimeResourcesViewList()
         {
-            return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_, particles_);
+            return std::forward_as_tuple(flux_, momentumTensor_, particleDensity_, chargeDensity_,
+                                         particles_);
         }