From 64a31a2e20a49bb08e595a7d5498c1c6a8d01f27 Mon Sep 17 00:00:00 2001
From: Nicolas Aunai <nicolas.aunai@lpp.polytechnique.fr>
Date: Tue, 21 Nov 2023 22:11:23 +0100
Subject: [PATCH] add massDensity

---
 src/core/data/ions/ions.hpp | 49 ++++++++++++++++++++++++++++++-------
 1 file changed, 40 insertions(+), 9 deletions(-)

diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp
index edd60dd23..d6e96627d 100644
--- a/src/core/data/ions/ions.hpp
+++ b/src/core/data/ions/ions.hpp
@@ -105,10 +105,37 @@ namespace core
                                std::begin(*rho_), std::plus<typename field_type::type>{});
             }
         }
+        void computeMassDensity()
+        {
+            massDensity_->zero();
+
+            for (auto const& pop : populations_)
+            {
+                // we sum over all nodes contiguously, including ghosts
+                // nodes. This is more efficient and easier to code as we don't
+                // have to account for the field dimensionality.
+
+                auto& popDensity = pop.density();
+                std::transform(
+                    std::begin(*rho_), std::end(*rho_), std::begin(popDensity), std::begin(*rho_),
+                    [&pop](auto const& n, auto const& pop_n) { return n + pop_n * pop.mass(); });
+            }
+        }
 
 
         void computeBulkVelocity()
         {
+            // at this point each population has its flux = density * velocity computed
+            // the bulk velocity is sum(pop_mass * pop_flux) / sum(pop_mass * pop_density)
+            // massDensity_ is computed by computeMassDensity() and stored in massDensity_
+            // rho_ is computed by computeDensity() already
+            // we can  use rho_ and not compute
+            // massDensity_ if all populations have the same mass since  sum(pop_mass * pop_flux) /
+            // sum(pop_mass * pop_density) would be equivalent to sum(pop_flux) / sum(pop_density)
+
+            computeMassDensity();
+            auto const& density = massDensity_;
+
             bulkVelocity_.zero();
             auto& vx = bulkVelocity_.getComponent(Component::X);
             auto& vy = bulkVelocity_.getComponent(Component::Y);
@@ -130,15 +157,12 @@ namespace core
                 std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz),
                                fluxSum);
             }
-            // auto denominator = [&](auto const& v, auto const& r) { return v / (pop.mass() * r);
-            // };
-
 
-            std::transform(std::begin(vx), std::end(vx), std::begin(*rho_), std::begin(vx),
+            std::transform(std::begin(vx), std::end(vx), std::begin(*massDensity_), std::begin(vx),
                            std::divides<typename field_type::type>{});
-            std::transform(std::begin(vy), std::end(vy), std::begin(*rho_), std::begin(vy),
+            std::transform(std::begin(vy), std::end(vy), std::begin(*massDensity_), std::begin(vy),
                            std::divides<typename field_type::type>{});
-            std::transform(std::begin(vz), std::end(vz), std::begin(*rho_), std::begin(vz),
+            std::transform(std::begin(vz), std::end(vz), std::begin(*massDensity_), std::begin(vz),
                            std::divides<typename field_type::type>{});
         }
 
@@ -152,6 +176,8 @@ namespace core
         NO_DISCARD auto end() const { return std::end(populations_); }
 
 
+        // in the following isUsable and isSettable the massDensity_ is not checked
+        // because it is for internal use only so no object will ever need to access it.
         NO_DISCARD bool isUsable() const
         {
             bool usable = rho_ != nullptr && bulkVelocity_.isUsable();
@@ -191,7 +217,8 @@ namespace core
 
         NO_DISCARD MomentProperties getFieldNamesAndQuantities() const
         {
-            return {{{densityName(), HybridQuantity::Scalar::rho}}};
+            return {{{densityName(), HybridQuantity::Scalar::rho},
+                     {"massDensity", HybridQuantity::Scalar::rho}}};
         }
 
 
@@ -202,6 +229,10 @@ namespace core
             {
                 rho_ = field;
             }
+            else if (bufferName == "massDensity")
+            {
+                massDensity_ = field;
+            }
             else
             {
                 throw std::runtime_error("Error - invalid density buffer name");
@@ -241,9 +272,9 @@ namespace core
 
     private:
         field_type* rho_{nullptr};
+        field_type* massDensity_{nullptr};
         vecfield_type bulkVelocity_;
-        std::vector<IonPopulation> populations_; // TODO we have to name this so they are unique
-                                                 // although only 1 Ions should exist.
+        std::vector<IonPopulation> populations_;
     };
 } // namespace core
 } // namespace PHARE