adding mass to bulk velocity (#783)

* account for population masses to bulk velocity * minor refac + compiler flag for uninitialized variable * have_sys_time --------- Co-authored-by: dekken <[email protected]>
PHAREHUB · Dec 2, 2023 · 61b1ade · 61b1ade
1 parent 501fb3b
commit 61b1ade
Show file tree

Hide file tree

Showing 21 changed files with 177 additions and 79 deletions.
diff --git a/.github/workflows/cmake_ubuntu.yml b/.github/workflows/cmake_ubuntu.yml
@@ -80,7 +80,7 @@ jobs:
         cmake $GITHUB_WORKSPACE -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON --fresh \
               -DCMAKE_C_COMPILER_LAUNCHER=ccache \
               -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \
-              -DlowResourceTests=ON \
+              -DlowResourceTests=ON -DdevMode=ON \
               -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 -O2"
 
     - name: Build

diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py
@@ -100,7 +100,10 @@ def add_int(path, val):
     def add_double(path, val):
         pp.add_double(path, float(val))
     def add_size_t(path, val):
-        pp.add_size_t(path, int(val))
+        casted = int(val)
+        if casted < 0:
+            raise RuntimeError("pyphare.__init__::add_size_t received negative value")
+        pp.add_size_t(path, casted)
     def add_vector_int(path, val):
         pp.add_vector_int(path, list(val))
 
@@ -192,8 +195,9 @@ def as_paths(rb):
     init_model = simulation.model
     modelDict  = init_model.model_dict
 
-    add_int("simulation/ions/nbrPopulations", init_model.nbr_populations())
-
+    if init_model.nbr_populations() < 0:
+        raise RuntimeError("Number of populations cannot be negative")
+    add_size_t("simulation/ions/nbrPopulations", init_model.nbr_populations())
 
     partinit = "particle_initializer"
     for pop_index, pop in enumerate(init_model.populations):

diff --git a/res/cmake/def.cmake b/res/cmake/def.cmake
@@ -1,6 +1,13 @@
 
 
-set (PHARE_FLAGS ${PHARE_FLAGS})
+# Per compiler CXXFLAGS
+set (PHARE_FLAGS ${PHARE_FLAGS} )
+if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+  set (PHARE_FLAGS ${PHARE_FLAGS} )
+else() # !Clang
+  set (PHARE_FLAGS ${PHARE_FLAGS} --param=min-pagesize=0 )
+endif() # clang
+
 set (PHARE_WERROR_FLAGS ${PHARE_FLAGS} ${PHARE_WERROR_FLAGS})
 set (PHARE_PYTHONPATH "${CMAKE_BINARY_DIR}:${CMAKE_SOURCE_DIR}/pyphare")
 set (PHARE_BASE_LIBS )
@@ -28,7 +35,7 @@ if(devMode) # -DdevMode=ON
   # Having quotes on strings here has lead to quotes being added to the compile string, so avoid.
 
   set (_Werr ${PHARE_WERROR_FLAGS} -Wall -Wextra -pedantic -Werror -Wno-unused-variable -Wno-unused-parameter)
-  set (_Werr ${_Werr} -Wdouble-promotion)
+  set (_Werr ${_Werr} -Wdouble-promotion -Wuninitialized )
 
   if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
     set (_Werr ${_Werr} -Wno-gnu-zero-variadic-macro-arguments)

diff --git a/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp b/src/amr/messengers/hybrid_hybrid_messenger_strategy.hpp
@@ -673,6 +673,7 @@ namespace amr
                                               NiOldUser_.name, fieldRefineOp_, fieldTimeOp_,
                                               info->modelIonDensity);
 
+
             velGhostsRefiners_.addTimeRefiners(info->ghostBulkVelocity, info->modelIonBulkVelocity,
                                                core::VecFieldNames{ViOld_}, fieldRefineOp_,
                                                fieldTimeOp_);

diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp
@@ -28,6 +28,7 @@ namespace core
     public:
         using field_type          = typename IonPopulation::field_type;
         using vecfield_type       = typename IonPopulation::vecfield_type;
+        using Float               = typename field_type::type;
         using particle_array_type = typename IonPopulation::particle_array_type;
         using ParticleInitializerFactoryT
             = ParticleInitializerFactory<particle_array_type, GridLayout>;
@@ -39,15 +40,13 @@ namespace core
 
         explicit Ions(PHARE::initializer::PHAREDict const& dict)
             : bulkVelocity_{"bulkVel", HybridQuantity::Vector::V}
-            , populations_{}
+            , populations_{generate(
+                  [&dict](auto ipop) { //
+                      return IonPopulation{dict["pop" + std::to_string(ipop)]};
+                  },
+                  dict["nbrPopulations"].template to<std::size_t>())}
+            , sameMasses_{allSameMass_()}
         {
-            auto nbrPop = dict["nbrPopulations"].template to<int>();
-
-            for (int ipop = 0; ipop < nbrPop; ++ipop)
-            {
-                auto const& pop = dict["pop" + std::to_string(ipop)];
-                populations_.emplace_back(pop);
-            }
         }
 
 
@@ -57,37 +56,29 @@ namespace core
         NO_DISCARD field_type const& density() const
         {
             if (isUsable())
-            {
                 return *rho_;
-            }
             else
-            {
                 throw std::runtime_error("Error - cannot access density data");
-            }
         }
 
 
 
         NO_DISCARD field_type& density()
         {
             if (isUsable())
-            {
                 return *rho_;
-            }
             else
-            {
                 throw std::runtime_error("Error - cannot access density data");
-            }
         }
 
 
 
-
         NO_DISCARD vecfield_type const& velocity() const { return bulkVelocity_; }
 
         NO_DISCARD vecfield_type& velocity() { return bulkVelocity_; }
 
         NO_DISCARD std::string densityName() const { return "rho"; }
+        NO_DISCARD std::string massDensityName() const { return "massDensity"; }
 
 
         void computeDensity()
@@ -102,43 +93,73 @@ namespace core
 
                 auto& popDensity = pop.density();
                 std::transform(std::begin(*rho_), std::end(*rho_), std::begin(popDensity),
-                               std::begin(*rho_), std::plus<typename field_type::type>{});
+                               std::begin(*rho_), std::plus<Float>{});
+            }
+        }
+        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(*massDensity_), std::end(*massDensity_), std::begin(popDensity),
+                    std::begin(*massDensity_),
+                    [&pop](auto const& n, auto const& pop_n) { return n + pop_n * pop.mass(); });
             }
         }
 
 
         void computeBulkVelocity()
         {
+            // the bulk velocity is sum(pop_mass * pop_flux) / sum(pop_mass * pop_density)
+            // if all populations have the same mass, this is equivalent to sum(pop_flux) /
+            // sum(pop_density) sum(pop_density) is rho_ and already known by the time we get here.
+            // sum(pop_mass * pop_flux) is massDensity_ and is computed by computeMassDensity() if
+            // needed
+            if (!sameMasses_)
+                computeMassDensity();
+            auto const& density = (sameMasses_) ? rho_ : massDensity_;
+
             bulkVelocity_.zero();
             auto& vx = bulkVelocity_.getComponent(Component::X);
             auto& vy = bulkVelocity_.getComponent(Component::Y);
             auto& vz = bulkVelocity_.getComponent(Component::Z);
 
             for (auto& pop : populations_)
             {
-                auto& flux = pop.flux();
-                auto& fx   = flux.getComponent(Component::X);
-                auto& fy   = flux.getComponent(Component::Y);
-                auto& fz   = flux.getComponent(Component::Z);
+                // account for mass only if populations have different masses
+                std::function<Float(Float, Float)> plus = std::plus<Float>{};
+                std::function<Float(Float, Float)> plusMass
+                    = [&pop](Float const& v, Float const& f) { return v + f * pop.mass(); };
+
+                auto const& flux    = pop.flux();
+                auto&& [fx, fy, fz] = flux();
+
 
                 std::transform(std::begin(vx), std::end(vx), std::begin(fx), std::begin(vx),
-                               std::plus<typename field_type::type>{});
+                               sameMasses_ ? plus : plusMass);
                 std::transform(std::begin(vy), std::end(vy), std::begin(fy), std::begin(vy),
-                               std::plus<typename field_type::type>{});
+                               sameMasses_ ? plus : plusMass);
                 std::transform(std::begin(vz), std::end(vz), std::begin(fz), std::begin(vz),
-                               std::plus<typename field_type::type>{});
+                               sameMasses_ ? plus : plusMass);
             }
 
-            std::transform(std::begin(vx), std::end(vx), std::begin(*rho_), std::begin(vx),
-                           std::divides<typename field_type::type>{});
-            std::transform(std::begin(vy), std::end(vy), std::begin(*rho_), std::begin(vy),
-                           std::divides<typename field_type::type>{});
-            std::transform(std::begin(vz), std::end(vz), std::begin(*rho_), std::begin(vz),
-                           std::divides<typename field_type::type>{});
+
+            std::transform(std::begin(vx), std::end(vx), std::begin(*density), std::begin(vx),
+                           std::divides<Float>{});
+            std::transform(std::begin(vy), std::end(vy), std::begin(*density), std::begin(vy),
+                           std::divides<Float>{});
+            std::transform(std::begin(vz), std::end(vz), std::begin(*density), std::begin(vz),
+                           std::divides<Float>{});
         }
 
 
-        // TODO 3347 compute ion bulk velocity
 
         NO_DISCARD auto begin() { return std::begin(populations_); }
         NO_DISCARD auto end() { return std::end(populations_); }
@@ -147,6 +168,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();
@@ -186,20 +209,25 @@ namespace core
 
         NO_DISCARD MomentProperties getFieldNamesAndQuantities() const
         {
-            return {{{densityName(), HybridQuantity::Scalar::rho}}};
+            return {{{densityName(), HybridQuantity::Scalar::rho},
+                     {massDensityName(), HybridQuantity::Scalar::rho}}};
         }
 
 
 
         void setBuffer(std::string const& bufferName, field_type* field)
         {
-            if (bufferName == "rho")
+            if (bufferName == densityName())
             {
                 rho_ = field;
             }
+            else if (bufferName == massDensityName())
+            {
+                massDensity_ = field;
+            }
             else
             {
-                throw std::runtime_error("Error - invalid density buffer name");
+                throw std::runtime_error("Error - invalid density buffer name : " + bufferName);
             }
         }
 
@@ -235,10 +263,25 @@ namespace core
 
 
     private:
+        bool allSameMass_() const
+        {
+            return all(populations_, [this](auto const& pop) { // arbitrary small diff
+                return float_equals(pop.mass(), populations_.front().mass(), /*abs_tol=*/1e-10);
+            });
+        }
+
+
+
+
         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_;
+
+        // note this is set at construction when all populations are added
+        // in the future if some populations are dynamically created during the simulation
+        // this should be updated accordingly
+        bool sameMasses_{false};
     };
 } // namespace core
 } // namespace PHARE

diff --git a/src/core/data/particles/particle_packer.hpp b/src/core/data/particles/particle_packer.hpp
@@ -28,7 +28,7 @@ class ParticlePacker
 
     static auto empty()
     {
-        Particle<dim> particle;
+        Particle<dim> particle{};
         return get(particle);
     }
 

diff --git a/src/core/utilities/types.hpp b/src/core/utilities/types.hpp
@@ -326,26 +326,34 @@ NO_DISCARD auto constexpr generate(F&& f, std::array<Type, Size> const& arr)
 auto constexpr static to_bool = [](auto const& v) { return bool{v}; };
 
 
-template<typename Container>
-NO_DISCARD auto all(Container const& container)
+template<typename Container, typename Fn = decltype(to_bool)>
+NO_DISCARD auto all(Container const& container, Fn fn = to_bool)
 {
-    return std::all_of(container.begin(), container.end(), to_bool);
+    return std::all_of(container.begin(), container.end(), fn);
 }
 
-template<typename Container>
-NO_DISCARD auto any(Container const& container)
+
+template<typename Container, typename Fn = decltype(to_bool)>
+NO_DISCARD auto any(Container const& container, Fn fn = to_bool)
 {
-    return std::any_of(container.begin(), container.end(), to_bool);
+    return std::any_of(container.begin(), container.end(), fn);
 }
 
-template<typename Container>
-NO_DISCARD auto none(Container const& container)
+template<typename Container, typename Fn = decltype(to_bool)>
+NO_DISCARD auto none(Container const& container, Fn fn = to_bool)
 {
-    return std::none_of(container.begin(), container.end(), to_bool);
+    return std::none_of(container.begin(), container.end(), fn);
 }
 
+auto inline float_equals(float const& a, float const& b, float diff = 1e-6)
+{
+    return std::abs(a - b) < diff;
+}
 
-
+auto inline float_equals(double const& a, double const& b, double diff = 1e-12)
+{
+    return std::abs(a - b) < diff;
+}
 
 } // namespace PHARE::core
 

diff --git a/src/initializer/python_data_provider.hpp b/src/initializer/python_data_provider.hpp
@@ -1,4 +1,3 @@
-
 // this file includes python which has conflicting #defines with SAMRAI
 //  so include it last
 
@@ -13,7 +12,7 @@
 DISABLE_WARNING(shadow, shadow-field-in-constructor-modified, 42)
 
 #undef HAVE_SYS_TIMES_H // included in python again, possibly with different value
-#undef HAVE_UNISTD_HPP
+#undef HAVE_UNISTD_H
 
 #include <pybind11/embed.h> // everything needed for embedding
 #include <pybind11/functional.h>

diff --git a/src/python3/CMakeLists.txt b/src/python3/CMakeLists.txt
@@ -4,7 +4,7 @@ project(phare_python3)
 
 pybind11_add_module(cpp cpp_simulator.cpp)
 target_link_libraries(cpp PUBLIC phare_simulator)
-target_compile_options(cpp PRIVATE ${PHARE_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) # pybind fails with Werror
+target_compile_options(cpp PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE})
 set_target_properties(cpp
     PROPERTIES
     LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs"
@@ -15,7 +15,7 @@ set_property(TARGET cpp PROPERTY INTERPROCEDURAL_OPTIMIZATION ${PHARE_INTERPROCE
 if (CMAKE_BUILD_TYPE STREQUAL "Debug")
   pybind11_add_module(cpp_dbg cpp_simulator.cpp)
   target_link_libraries(cpp_dbg PUBLIC phare_simulator)
-  target_compile_options(cpp_dbg PRIVATE ${PHARE_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg)
+  target_compile_options(cpp_dbg PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE} -DPHARE_DIAG_DOUBLES=1 -DPHARE_CPP_MOD_NAME=cpp_dbg)
   set_target_properties(cpp_dbg
       PROPERTIES
       LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs"
@@ -26,7 +26,7 @@ endif (CMAKE_BUILD_TYPE STREQUAL "Debug")
 
 pybind11_add_module(cpp_etc cpp_etc.cpp)
 target_link_libraries(cpp_etc PUBLIC phare_simulator)
-target_compile_options(cpp_etc PRIVATE ${PHARE_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) # pybind fails with Werror
+target_compile_options(cpp_etc PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE})
 set_target_properties(cpp_etc
     PROPERTIES
     LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/pybindlibs"

diff --git a/src/python3/cpp_simulator.hpp b/src/python3/cpp_simulator.hpp
@@ -13,13 +13,13 @@
 #include "phare/phare.hpp"
 #include "simulator/simulator.hpp"
 
+#include "python3/pybind_def.hpp"
 #include "pybind11/stl.h"
 #include "pybind11/numpy.h"
 #include "pybind11/chrono.h"
 #include "pybind11/complex.h"
 #include "pybind11/functional.h"
 
-#include "python3/pybind_def.hpp"
 #include "python3/particles.hpp"
 #include "python3/patch_data.hpp"
 #include "python3/patch_level.hpp"