add mass density diag and refac compute_hier_from

pyphare/pyphare/pharein/diagnostics.py

@@ -205,7 +205,13 @@ def population_in_model(population):
 
 class FluidDiagnostics_(Diagnostics):
 
-    fluid_quantities = ["density", "flux", "bulkVelocity", "momentum_tensor"]
+    fluid_quantities = [
+        "density",
+        "mass_density",
+        "flux",
+        "bulkVelocity",
+        "momentum_tensor",
+    ]
     type = "fluid"
 
     def __init__(self, **kwargs):
@@ -268,7 +274,7 @@ class FluidDiagnostics:
     def __init__(self, **kwargs):
         if kwargs["quantity"] == "pressure_tensor":
             if for_total_ions(**kwargs):
-                needed_quantities = ["density", "bulkVelocity", "momentum_tensor"]
+                needed_quantities = ["mass_density", "bulkVelocity", "momentum_tensor"]
             else:
                 needed_quantities = ["density", "flux", "momentum_tensor"]
 

pyphare/pyphare/pharesee/hierarchy.py

@@ -1169,6 +1169,7 @@ def is_root_lvl(patch_level):
     "momentum_tensor_yz": "Myz",
     "momentum_tensor_zz": "Mzz",
     "density": "rho",
+    "mass_density": "rho",
     "tags": "tags",
 }
 
@@ -1256,7 +1257,60 @@ def compute_hier_from(h, compute):
 
 
 def are_compatible_hierarchies(hierarchies):
-    pass
+    return True
+
+
+def extract_patchdatas(hierarchies, ilvl, ipatch):
+    """
+    returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl
+    """
+    patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies]
+    patch_datas = {
+        pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items())
+    }
+    return patch_datas
+
+
+def new_patchdatas_from(compute, patchdatas, layout):
+    new_patch_datas = {}
+    datas = compute(patchdatas)
+    for data in datas:
+        pd = FieldData(layout, data["name"], data["data"], centering=data["centering"])
+        new_patch_datas[data["name"]] = pd
+    return new_patch_datas
+
+
+def new_patches_from(compute, hierarchies, ilvl):
+    reference_hier = hierarchies[0]
+    new_patches = []
+    patch_nbr = len(reference_hier.patch_levels[ilvl].patches)
+    for ip in range(patch_nbr):
+        current_patch = reference_hier.patch_levels[ilvl].patches[ip]
+        layout = current_patch.layout
+        patch_datas = extract_patchdatas(hierarchies, ilvl, ip)
+        new_patch_datas = new_patchdatas_from(compute, patch_datas, layout)
+        new_patches.append(Patch(new_patch_datas, current_patch.id))
+    return new_patches
+
+
+def compute_hier_from__(compute, hierarchies):
+    """return a new hierarchy using the callback 'compute' on all patchdatas
+    of the given hierarchies
+    """
+    if not are_compatible_hierarchies(hierarchies):
+        raise RuntimeError("hierarchies are not compatible")
+
+    reference_hier = hierarchies[0]
+    domain_box = reference_hier.domain_box
+    patch_levels = {}
+    for ilvl in range(reference_hier.levelNbr()):
+        patch_levels[ilvl] = PatchLevel(
+            ilvl, new_patches_from(compute, hierarchies, ilvl)
+        )
+
+    assert len(reference_hier.time_hier) == 1  # only single time hierarchies now
+    t = list(reference_hier.time_hier.keys())[0]
+    return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t)
 
 
 def compute_hier_from_(compute, hierarchies):

src/core/data/ions/ions.hpp

@@ -73,6 +73,23 @@ namespace core
                 throw std::runtime_error("Error - cannot access density data");
         }
 
+        NO_DISCARD field_type const& massDensity() const
+        {
+            if (isUsable())
+                return sameMasses_ ? *rho_ : *massDensity_;
+            else
+                throw std::runtime_error("Error - cannot access density data");
+        }
+
+
+
+        NO_DISCARD field_type& massDensity()
+        {
+            if (isUsable())
+                return sameMasses_ ? *rho_ : *massDensity_;
+            else
+                throw std::runtime_error("Error - cannot access density data");
+        }
 
 
         NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; }
@@ -195,6 +212,10 @@ namespace core
         {
             bool usable
                 = rho_ != nullptr and bulkVelocity_.isUsable() and momentumTensor_.isUsable();
+
+            // if all populations have the same mass, we don't need the massDensity_
+            usable &= (sameMasses_) ? true : massDensity_ != nullptr;
+
             for (auto const& pop : populations_)
             {
                 usable = usable and pop.isUsable();
@@ -208,6 +229,10 @@ namespace core
         {
             bool settable
                 = rho_ == nullptr and bulkVelocity_.isSettable() and momentumTensor_.isSettable();
+
+            // if all populations have the same mass, we don't need the massDensity_
+            settable &= (sameMasses_) ? true : massDensity_ == nullptr;
+
             for (auto const& pop : populations_)
             {
                 settable = settable and pop.isSettable();
@@ -228,12 +253,15 @@ namespace core
             typename HybridQuantity::Scalar qty;
         };
 
-        using MomentProperties = std::array<MomentsProperty, 2>;
+        using MomentProperties = std::vector<MomentsProperty>;
 
         NO_DISCARD MomentProperties getFieldNamesAndQuantities() const
         {
-            return {{{densityName(), HybridQuantity::Scalar::rho},
-                     {massDensityName(), HybridQuantity::Scalar::rho}}};
+            if (sameMasses_)
+                return {{{densityName(), HybridQuantity::Scalar::rho}}};
+            else
+                return {{{densityName(), HybridQuantity::Scalar::rho},
+                         {massDensityName(), HybridQuantity::Scalar::rho}}};
         }
 
 
@@ -246,6 +274,7 @@ namespace core
             }
             else if (bufferName == massDensityName())
             {
+                assert(sameMasses_ == false);
                 massDensity_ = field;
             }
             else
@@ -284,6 +313,8 @@ namespace core
             return ss.str();
         }
 
+        NO_DISCARD bool sameMasses() const { return sameMasses_; }
+
 
     private:
         bool allSameMass_() const

src/diagnostic/detail/types/fluid.hpp

@@ -100,9 +100,9 @@ void FluidDiagnosticWriter<H5Writer>::compute(DiagnosticProperties& diagnostic)
                       // and the former levelGhostPartsNew has been moved to levelGhostPartsOld
                       auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesOld());
 
-                      interpolator(domainParts, pop_momentum_tensor, layout);
-                      interpolator(patchGhostParts, pop_momentum_tensor, layout);
-                      interpolator(levelGhostParts, pop_momentum_tensor, layout);
+                      interpolator(domainParts, pop_momentum_tensor, layout, pop.mass());
+                      interpolator(patchGhostParts, pop_momentum_tensor, layout, pop.mass());
+                      interpolator(levelGhostParts, pop_momentum_tensor, layout, pop.mass());
                   }
                   ions.computeFullMomentumTensor();
               };
@@ -126,9 +126,9 @@ void FluidDiagnosticWriter<H5Writer>::compute(DiagnosticProperties& diagnostic)
                       // and the former levelGhostPartsNew has been moved to levelGhostPartsOld
                       auto levelGhostParts = core::makeIndexRange(pop.levelGhostParticlesOld());
 
-                      interpolator(domainParts, pop_momentum_tensor, layout);
-                      interpolator(patchGhostParts, pop_momentum_tensor, layout);
-                      interpolator(levelGhostParts, pop_momentum_tensor, layout);
+                      interpolator(domainParts, pop_momentum_tensor, layout, pop.mass());
+                      interpolator(patchGhostParts, pop_momentum_tensor, layout, pop.mass());
+                      interpolator(levelGhostParts, pop_momentum_tensor, layout, pop.mass());
                   };
             if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
             {
@@ -150,7 +150,8 @@ void FluidDiagnosticWriter<H5Writer>::createFiles(DiagnosticProperties& diagnost
     }
 
     std::string tree{"/ions/"};
-    checkCreateFileFor_(diagnostic, fileData_, tree, "density", "bulkVelocity", "momentum_tensor");
+    checkCreateFileFor_(diagnostic, fileData_, tree, "density", "mass_density", "bulkVelocity",
+                        "momentum_tensor");
 }
 
 
@@ -207,6 +208,8 @@ void FluidDiagnosticWriter<H5Writer>::getDataSetInfo(DiagnosticProperties& diagn
     std::string tree{"/ions/"};
     if (isActiveDiag(diagnostic, tree, "density"))
         infoDS(ions.density(), "density", patchAttributes[lvlPatchID]["ion"]);
+    if (isActiveDiag(diagnostic, tree, "mass_density"))
+        infoDS(ions.massDensity(), "mass_density", patchAttributes[lvlPatchID]["ion"]);
     if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
         infoVF(ions.velocity(), "bulkVelocity", patchAttributes[lvlPatchID]["ion"]);
     if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -274,6 +277,8 @@ void FluidDiagnosticWriter<H5Writer>::initDataSets(
         std::string tree{"/ions/"};
         if (isActiveDiag(diagnostic, tree, "density"))
             initDS(path, attr["ion"], "density", null);
+        if (isActiveDiag(diagnostic, tree, "mass_density"))
+            initDS(path, attr["ion"], "mass_density", null);
         if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
             initVF(path, attr["ion"], "bulkVelocity", null);
         if (isActiveDiag(diagnostic, tree, "momentum_tensor"))
@@ -310,9 +315,10 @@ void FluidDiagnosticWriter<H5Writer>::write(DiagnosticProperties& diagnostic)
     }
 
     std::string tree{"/ions/"};
-    auto& density = ions.density();
     if (isActiveDiag(diagnostic, tree, "density"))
-        writeDS(path + "density", density);
+        writeDS(path + "density", ions.density());
+    if (isActiveDiag(diagnostic, tree, "mass_density"))
+        writeDS(path + "mass_density", ions.massDensity());
     if (isActiveDiag(diagnostic, tree, "bulkVelocity"))
         writeTF(path + "bulkVelocity", ions.velocity());
     if (isActiveDiag(diagnostic, tree, "momentum_tensor"))

tests/core/numerics/ion_updater/test_updater.cpp

@@ -466,7 +466,8 @@ struct IonsBuffers
     void setBuffers(Ions& ions)
     {
         ions.setBuffer("rho", &ionDensity);
-        ions.setBuffer("massDensity", &ionMassDensity);
+        if (!ions.sameMasses())
+            ions.setBuffer("massDensity", &ionMassDensity);
         auto& v = ions.velocity();
         v.setBuffer("bulkVel_x", &Vx);
         v.setBuffer("bulkVel_y", &Vy);