introduce charge, mass and particle densities

.gdbinit

@@ -0,0 +1,3 @@
+break __assert_fail
+break abort
+catch throw

.gitignore

@@ -22,3 +22,4 @@ perf.*
 .vscode
 .phare*
 PHARE_REPORT.zip
+.gdbinit

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,7 @@ 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:
+        return hierarchy_from_func(func, hier, **kwargs)
+
     raise ValueError("can't make hierarchy")

pyphare/pyphare/pharesee/hierarchy/fromfunc.py

@@ -0,0 +1,97 @@
+from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from
+
+import numpy as np
+
+
+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
+
+    assert len(masses) == len(densities)
+    funcs  = np.zeros((x.size, len(masses)))
+
+    for i, (mass, density) in enumerate(zip(masses, densities)):
+        funcs[:,i] = mass*density(x)
+
+    return funcs.sum(axis=1)
+
+
+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)
+
+    funcs  = np.zeros((x.size, len(charges)))
+
+    for i, (charge, density) in enumerate(zip(charges, densities)):
+        funcs[:,i] = charge*density(x)
+
+    return funcs.sum(axis=1)
+
+
+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 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
+
+    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)):
+        funcs[:,:,i] = mass*density(xv, yv)
+
+    return funcs.sum(axis=2)
+
+
+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
+
+    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)):
+        funcs[:,:,i] = charge*density(xv, yv)
+
+    return funcs.sum(axis=2)
+
+
+def hierarchy_from_func2d(func, hier, **kwargs):
+    assert hier.ndim == 2
+
+    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):
+    if hier.ndim == 1:
+        return hierarchy_from_func1d(func, hier, **kwargs)
+    if hier.ndim == 2:
+        return hierarchy_from_func2d(func, hier, **kwargs)
+

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -360,6 +360,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]
         ):
@@ -370,35 +372,47 @@ 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)))
+
+        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
-                )
-            )
+            value = collection.get("value", np.nan)
+            for box in collection["boxes"]:
+                i0, j0 = box.lower
+                i1, j1 = box.upper
+                ij[i0 : i1 + 1, j0 : j1 + 1] = value
+            if "coords" in collection:
+                for coords in collection["coords"]:
+                    ij[coords] = collection["value"]
+
+        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
 
@@ -433,7 +447,7 @@ def plot1d(self, **kwargs):
                     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)
@@ -544,9 +558,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",
+                    "value": 1,
                 },
                 {
                     "boxes": [p.box for p in hierarchy.level(ilvl).patches],
-                    "facecolor": "grey",
+                    "value": 2,
                 },
             ],
         )
@@ -390,11 +390,11 @@ def test_level_ghost_boxes(self, refinement_boxes, expected):
                 collections=[
                     {
                         "boxes": ghost_area_box_list,
-                        "facecolor": "yellow",
+                        "value": 1,
                     },
                     {
                         "boxes": [p.box for p in hierarchy.level(ilvl).patches],
-                        "facecolor": "grey",
+                        "value": 2,
                     },
                 ],
                 title="".join(
@@ -502,7 +502,7 @@ def test_patch_periodicity_copy(self, refinement_boxes, expected):
                 collections=[
                     {
                         "boxes": [p.box for p in periodic_list],
-                        "facecolor": "grey",
+                        "value": 2,
                     },
                 ],
             )

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 @@
     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, _]