Charge

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/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",
 }
 

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",
 }
 
@@ -116,7 +116,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):
@@ -153,7 +153,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_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

@@ -34,7 +34,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"]}
         {
@@ -45,21 +46,18 @@ 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 particles_.isUsable() && rho_.isUsable() && flux_.isUsable()
-                   && momentumTensor_.isUsable();
+            return particles_.isUsable() && particleDensity_.isUsable() && chargeDensity_.isUsable()
+                   && flux_.isUsable() && momentumTensor_.isUsable();
         }
 
-
         NO_DISCARD bool isSettable() const
         {
-            return particles_.isSettable() && rho_.isSettable() && flux_.isSettable()
+            return particles_.isSettable() && particleDensity_.isSettable()
+                   && chargeDensity_.isSettable() && flux_.isSettable()
                    && momentumTensor_.isSettable();
         }
 
@@ -84,9 +82,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_; }
@@ -105,7 +105,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_);
         }
 
 
@@ -129,7 +130,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_;
     };