wip

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

@@ -363,7 +363,7 @@ 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)))
 
         for collection in collections:
             facecolor = collection.get("facecolor", "none")
@@ -383,15 +383,23 @@ def plot_2d_patches(self, ilvl, collections, **kwargs):
             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)
+        major_ticks = np.arange(mi - 5, ma + 5 + 5, 1)
         ax.set_xticks(major_ticks)
         ax.set_yticks(major_ticks)
 
-        minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1)
+        minor_ticks = np.arange(mi - 4.5, ma + 5 + 5, 1)
         ax.set_xticks(minor_ticks, minor=True)
         ax.set_yticks(minor_ticks, minor=True)
 
+        # Align the minor tick label
+        for label in ax.get_xticklabels(minor=True):
+            label.set_horizontalalignment("center")
+
         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/diagnostic/detail/types/fluid.hpp

@@ -199,7 +199,9 @@ void FluidDiagnosticWriter<H5Writer>::getDataSetInfo(DiagnosticProperties& diagn
         std::string tree{"/ions/pop/" + pop.name() + "/"};
         auto& popAttr = patchAttributes[lvlPatchID]["fluid_" + pop.name()];
         if (isActiveDiag(diagnostic, tree, "density"))
-            infoDS(pop.chargeDensity(), "density", popAttr);
+            infoDS(pop.particleDensity(), "density", popAttr);
+        if (isActiveDiag(diagnostic, tree, "charge_density"))
+            infoDS(pop.chargeDensity(), "charge_density", popAttr);
         if (isActiveDiag(diagnostic, tree, "flux"))
             infoVF(pop.flux(), "flux", popAttr);
         if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -269,6 +271,8 @@ void FluidDiagnosticWriter<H5Writer>::initDataSets(
             std::string popPath(path + "pop/" + pop.name() + "/");
             if (isActiveDiag(diagnostic, tree, "density"))
                 initDS(path, attr[popId], "density", null);
+            if (isActiveDiag(diagnostic, tree, "charge_density"))
+                initDS(path, attr[popId], "charge_density", null);
             if (isActiveDiag(diagnostic, tree, "flux"))
                 initVF(path, attr[popId], "flux", null);
             if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -308,7 +312,9 @@ void FluidDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
     {
         std::string tree{"/ions/pop/" + pop.name() + "/"};
         if (isActiveDiag(diagnostic, tree, "density"))
-            writeDS(path + "density", pop.chargeDensity());
+            writeDS(path + "density", pop.particleDensity());
+        if (isActiveDiag(diagnostic, tree, "charge_density"))
+            writeDS(path + "charge_density", pop.chargeDensity());
         if (isActiveDiag(diagnostic, tree, "flux"))
             writeTF(path + "flux", pop.flux());
         if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -345,6 +351,7 @@ void FluidDiagnosticWriter<H5Writer>::writeAttributes(
     for (auto& pop : h5Writer.modelView().getIons())
     {
         std::string tree = "/ions/pop/" + pop.name() + "/";
+        checkWrite(tree, "density", pop);
         checkWrite(tree, "charge_density", pop);
         checkWrite(tree, "flux", pop);
         checkWrite(tree, "momentum_tensor", pop);

tests/diagnostic/__init__.py

@@ -28,15 +28,20 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None):
     #     write_timestamps=timestamps
     # )
 
-    for quantity in ["charge_density", "bulkVelocity", "pressure_tensor"]:
+    for quantity in [
+        "charge_density",
+        "mass_density",
+        "bulkVelocity",
+        "pressure_tensor",
+    ]:
         ph.FluidDiagnostics(
             quantity=quantity,
             write_timestamps=timestamps,
             flush_every=flush_every,
         )
 
     for pop in pops:
-        for quantity in ["density", "flux", "pressure_tensor"]:
+        for quantity in ["density", "charge_density", "flux", "pressure_tensor"]:
             ph.FluidDiagnostics(
                 quantity=quantity,
                 write_timestamps=timestamps,

tests/diagnostic/test_diagnostics.hpp

@@ -226,7 +226,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
     using GridLayout                           = typename Simulator::PHARETypes::GridLayout_t;
     constexpr auto dimension                   = Simulator::dimension;
     constexpr std::size_t expectedPopNbr       = 2;
-    constexpr std::size_t expectedPopAttrFiles = 5;
+    constexpr std::size_t expectedPopAttrFiles = 6;
 
     std::string const ionsPopPath = "/ions/pop/";
 
@@ -236,7 +236,8 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
     auto nbrPop = dict["simulation"]["ions"]["nbrPopulations"].template to<std::size_t>();
     EXPECT_EQ(nbrPop, expectedPopNbr);
 
-    std::vector<std::string> h5FileTypes{"/EM_B", "/EM_E", "/ions/charge_density", "/ions/bulkVelocity"};
+    std::vector<std::string> h5FileTypes{"/EM_B", "/EM_E", "/ions/charge_density",
+                                         "/ions/mass_density", "/ions/bulkVelocity"};
 
     for (std::size_t i = 0; i < nbrPop; ++i)
     {
@@ -247,6 +248,7 @@ void validateAttributes(Simulator& sim, Hi5Diagnostic& hi5)
         h5FileTypes.emplace_back(ionsPopPath + popName + "/levelGhost");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/patchGhost");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/density");
+        h5FileTypes.emplace_back(ionsPopPath + popName + "/charge_density");
         h5FileTypes.emplace_back(ionsPopPath + popName + "/flux");
     }
 

tests/simulator/test_advance.py

@@ -47,7 +47,7 @@ def getHierarchy(
         largest_patch_size=20,
         cells=120,
         time_step=0.001,
-        model_init={},
+        model_init={"seed": 2},
         dl=0.2,
         extra_diag_options={},
         time_step_nbr=1,
@@ -289,10 +289,25 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time):
                         # seems correct considering ghosts are filled with schedules
                         # involving linear/spatial interpolations and so on where
                         # rounding errors may occur.... setting atol to 5.5e-15
-                        assert_fp_any_all_close(slice1, slice2, atol=5.5e-15, rtol=0)
+                        assert_fp_any_all_close(slice1, slice2, atol=2.5e-14, rtol=0)
                         checks += 1
                     except AssertionError as e:
                         print("AssertionError", pd1.name, e)
+                        print(
+                            np.where(~np.isclose(slice1, slice2, atol=2.5e-14, rtol=0))
+                        )
+
+                        fig = datahier.plot_2d_patches(
+                            ilvl,
+                            collections=[
+                                {"boxes": [pd1.box], "facecolor": "red"},
+                                {"boxes": [pd2.box], "facecolor": "blue"},
+                            ],
+                            xlim=(-5, 50),
+                            ylim=(-5, 50),
+                        )
+                        fig.savefig(f"pd1.png")
+                        print(f"ilvl {ilvl}")
                         print(pd1.box, pd2.box)
                         print(pd1.x.mean())
                         print(pd1.y.mean())