try float32

res/cmake/dep/highfive.cmake

@@ -4,9 +4,9 @@ set (PHARE_HAS_HIGHFIVE "0")
 if(HighFive)
 
   set (HIGHFIVE_SRC ${CMAKE_CURRENT_SOURCE_DIR}/subprojects/highfive)
-  set (HIGHFIVE_VERSION master)
+  set (HIGHFIVE_VERSION main)
 
-  phare_github_get_or_update(HighFive ${HIGHFIVE_SRC} BlueBrain/HighFive ${HIGHFIVE_VERSION})
+  phare_github_get_or_update(HighFive ${HIGHFIVE_SRC} highfive-devs/highfive ${HIGHFIVE_VERSION})
 
   include_directories(
     ${HIGHFIVE_SRC}/include

src/core/CMakeLists.txt

@@ -51,7 +51,6 @@ set( SOURCES_INC
    )
 
 set( SOURCES_CPP
-     data/ions/particle_initializers/maxwellian_particle_initializer.cpp
      utilities/index/index.cpp
      utilities/mpi_utils.cpp
     )

src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp

@@ -1,32 +1,24 @@
 #ifndef PHARE_FLUID_PARTICLE_INITIALIZER_HPP
 #define PHARE_FLUID_PARTICLE_INITIALIZER_HPP
 
-#include <memory>
-#include <random>
-#include <cassert>
-#include <functional>
 
-#include "core/data/grid/gridlayoutdefs.hpp"
-#include "core/hybrid/hybrid_quantities.hpp"
+#include "core/def.hpp"
 #include "core/utilities/types.hpp"
-#include "core/data/ions/particle_initializers/particle_initializer.hpp"
-#include "core/data/particles/particle.hpp"
 #include "initializer/data_provider.hpp"
 #include "core/utilities/point/point.hpp"
-#include "core/def.hpp"
-
-
-namespace PHARE::core
-{
-void maxwellianVelocity(std::array<double, 3> const& V, std::array<double, 3> const& Vth,
-                        std::mt19937_64& generator, std::array<double, 3>& partVelocity);
+#include "core/data/particles/particle.hpp"
+#include "core/data/grid/gridlayoutdefs.hpp"
+#include "core/data/ions/particle_initializers/particle_initializer.hpp"
 
+#include "maxwellian_particle_initializer_base.hpp"
 
-std::array<double, 3> basisTransform(std::array<std::array<double, 3>, 3> const& basis,
-                                     std::array<double, 3> const& vec);
+#include <memory>
+#include <random>
+#include <cassert>
 
-void localMagneticBasis(std::array<double, 3> B, std::array<std::array<double, 3>, 3>& basis);
 
+namespace PHARE::core
+{
 
 /** @brief a MaxwellianParticleInitializer is a ParticleInitializer that loads particles from a
  * local Maxwellian distribution given density, bulk velocity and thermal velocity profiles.
@@ -148,7 +140,7 @@ void MaxwellianParticleInitializer<ParticleArray, GridLayout>::loadParticles(
     };
 
 
-    auto deltas = [](auto& pos, auto& gen) -> std::array<double, dimension> {
+    auto deltas = [](auto& pos, auto& gen) -> std::array<floater_t<0>, dimension> {
         if constexpr (dimension == 1)
             return {pos(gen)};
         if constexpr (dimension == 2)
@@ -176,7 +168,7 @@ void MaxwellianParticleInitializer<ParticleArray, GridLayout>::loadParticles(
 
     auto const [n, V, Vth] = fns();
     auto randGen           = getRNG(rngSeed_);
-    ParticleDeltaDistribution<double> deltaDistrib;
+    ParticleDeltaDistribution<floater_t<0>> deltaDistrib;
 
     for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); ++flatCellIdx)
     {
@@ -186,20 +178,20 @@ void MaxwellianParticleInitializer<ParticleArray, GridLayout>::loadParticles(
         auto const cellWeight   = n[flatCellIdx] / nbrParticlePerCell_;
         auto const AMRCellIndex = layout.localToAMR(point(flatCellIdx, ndCellIndices));
         auto const iCell        = AMRCellIndex.template toArray<int>();
-        std::array<double, 3> particleVelocity;
-        std::array<std::array<double, 3>, 3> basis;
+        std::array<floater_t<1>, 3> particleVelocity;
+        std::array<std::array<double, 3>, 3> basis{};
 
         if (basis_ == Basis::Magnetic)
         {
             auto const B = fns.B();
-            localMagneticBasis({B[0][flatCellIdx], B[1][flatCellIdx], B[2][flatCellIdx]}, basis);
+            localMagneticBasis(basis, B[0][flatCellIdx], B[1][flatCellIdx], B[2][flatCellIdx]);
         }
 
         for (std::uint32_t ipart = 0; ipart < nbrParticlePerCell_; ++ipart)
         {
-            maxwellianVelocity({V[0][flatCellIdx], V[1][flatCellIdx], V[2][flatCellIdx]},
-                               {Vth[0][flatCellIdx], Vth[1][flatCellIdx], Vth[2][flatCellIdx]}, //
-                               randGen, particleVelocity);
+            maxwellianVelocity(particleVelocity, randGen, V[0][flatCellIdx], V[1][flatCellIdx],
+                               V[2][flatCellIdx], Vth[0][flatCellIdx], Vth[1][flatCellIdx],
+                               Vth[2][flatCellIdx]);
 
             if (basis_ == Basis::Magnetic)
                 particleVelocity = basisTransform(basis, particleVelocity);

src/core/data/ions/particle_initializers/maxwellian_particle_initializer.cpp → src/core/data/ions/particle_initializers/maxwellian_particle_initializer_base.hpp

@@ -1,20 +1,28 @@
+#ifndef PHARE_MAXWELLIAN_PARTICLE_INITIALIZER_BASE_HPP
+#define PHARE_MAXWELLIAN_PARTICLE_INITIALIZER_BASE_HPP
+
 #include <array>
+#include <tuple>
 #include <random>
 
-#include "core/utilities/types.hpp"
 #include "core/def.hpp"
+#include "core/utilities/types.hpp"
 
 
 namespace PHARE
 {
 namespace core
 {
-    void maxwellianVelocity(std::array<double, 3> const& V, std::array<double, 3> const& Vth,
-                            std::mt19937_64& generator, std::array<double, 3>& partVelocity)
+
+    template<typename T, typename... Args>
+    void maxwellianVelocity(std::array<T, 3>& partVelocity, std::mt19937_64& generator,
+                            Args const... args)
     {
-        std::normal_distribution<> maxwellX(V[0], Vth[0]);
-        std::normal_distribution<> maxwellY(V[1], Vth[1]);
-        std::normal_distribution<> maxwellZ(V[2], Vth[2]);
+        auto const& [V0, V1, V2, Vth0, Vth1, Vth2] = std::forward_as_tuple(args...);
+
+        std::normal_distribution<> maxwellX(V0, Vth0);
+        std::normal_distribution<> maxwellY(V1, Vth1);
+        std::normal_distribution<> maxwellZ(V2, Vth2);
 
         partVelocity[0] = maxwellX(generator);
         partVelocity[1] = maxwellY(generator);
@@ -23,12 +31,11 @@ namespace core
 
 
 
-
-    NO_DISCARD std::array<double, 3>
-    basisTransform(std::array<std::array<double, 3>, 3> const& basis,
-                   std::array<double, 3> const& vec)
+    template<typename T>
+    NO_DISCARD std::array<T, 3> basisTransform(std::array<std::array<double, 3>, 3> const& basis,
+                                               std::array<T, 3> const& vec)
     {
-        std::array<double, 3> newVec;
+        std::array<T, 3> newVec;
 
         for (std::uint32_t comp = 0; comp < 3; comp++)
         {
@@ -40,9 +47,11 @@ namespace core
     }
 
 
-
-    void localMagneticBasis(std::array<double, 3> B, std::array<std::array<double, 3>, 3>& basis)
+    template<typename T0, typename... Bargs>
+    void localMagneticBasis(std::array<std::array<T0, 3>, 3>& basis, Bargs const... bargs)
     {
+        std::array const B{bargs...};
+
         auto b2 = norm(B);
 
         if (b2 < 1e-8)
@@ -85,3 +94,5 @@ namespace core
 
 } // namespace core
 } // namespace PHARE
+
+#endif /*PHARE_MAXWELLIAN_PARTICLE_INITIALIZER_BASE_HPP*/

src/core/data/ndarray/ndarray_vector.hpp

@@ -2,16 +2,16 @@
 #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP
 
 #include "core/def.hpp"
-#include <stdexcept>
+
+#include "core/utilities/types.hpp"
+
+// #include <iostream>
 #include <array>
-#include <cstdint>
-#include <vector>
 #include <tuple>
+#include <vector>
+#include <cstdint>
 #include <numeric>
-#include <iostream>
-
-
-#include "core/utilities/types.hpp"
+#include <stdexcept>
 
 
 namespace PHARE::core

src/core/data/particles/particle.hpp

@@ -20,12 +20,15 @@ namespace PHARE::core
 template<typename T = float>
 struct ParticleDeltaDistribution
 {
+    T constexpr static one = 1;
+    T constexpr static zro = 0;
+
     template<typename Generator>
     NO_DISCARD T operator()(Generator& generator)
     {
         return dist(generator);
     }
-    std::uniform_real_distribution<T> dist{0, 1. - std::numeric_limits<T>::epsilon()};
+    std::uniform_real_distribution<T> dist{zro, one - std::numeric_limits<T>::epsilon()};
 };
 
 
@@ -42,8 +45,8 @@ struct Particle
     static_assert(dim > 0 and dim < 4, "Only dimensions 1,2,3 are supported.");
     static const size_t dimension = dim;
 
-    Particle(double a_weight, double a_charge, std::array<int, dim> cell,
-             std::array<double, dim> a_delta, std::array<double, 3> a_v)
+    Particle(floater_t<3> a_weight, floater_t<2> a_charge, std::array<int, dim> cell,
+             std::array<floater_t<0>, dim> a_delta, std::array<floater_t<1>, 3> a_v)
         : weight{a_weight}
         , charge{a_charge}
         , iCell{cell}
@@ -54,12 +57,12 @@ struct Particle
 
     Particle() = default;
 
-    double weight = 0;
-    double charge = 0;
+    floater_t<3> weight = 0;
+    floater_t<2> charge = 0;
 
-    std::array<int, dim> iCell    = ConstArray<int, dim>();
-    std::array<double, dim> delta = ConstArray<double, dim>();
-    std::array<double, 3> v       = ConstArray<double, 3>();
+    std::array<int, dim> iCell          = ConstArray<int, dim>();
+    std::array<floater_t<0>, dim> delta = ConstArray<floater_t<0>, dim>();
+    std::array<floater_t<1>, 3> v       = ConstArray<floater_t<1>, 3>();
 
     NO_DISCARD bool operator==(Particle<dim> const& that) const
     {
@@ -104,11 +107,11 @@ struct ParticleView
     static_assert(dim > 0 and dim < 4, "Only dimensions 1,2,3 are supported.");
     static constexpr std::size_t dimension = dim;
 
-    double& weight;
-    double& charge;
+    floater_t<3>& weight;
+    floater_t<2>& charge;
     std::array<int, dim>& iCell;
-    std::array<double, dim>& delta;
-    std::array<double, 3>& v;
+    std::array<floater_t<0>, dim>& delta;
+    std::array<floater_t<1>, 3>& v;
 };
 
 
@@ -121,9 +124,9 @@ inline constexpr auto is_phare_particle_type
 
 template<std::size_t dim, template<std::size_t> typename ParticleA,
          template<std::size_t> typename ParticleB>
-NO_DISCARD typename std::enable_if_t<
-    is_phare_particle_type<dim, ParticleA<dim>> and is_phare_particle_type<dim, ParticleB<dim>>,
-    bool>
+NO_DISCARD typename std::enable_if_t<is_phare_particle_type<dim, ParticleA<dim>>
+                                         and is_phare_particle_type<dim, ParticleB<dim>>,
+                                     bool>
 operator==(ParticleA<dim> const& particleA, ParticleB<dim> const& particleB)
 {
     return particleA.weight == particleB.weight and //

src/core/data/particles/particle_array.hpp

@@ -2,18 +2,18 @@
 #define PHARE_CORE_DATA_PARTICLES_PARTICLE_ARRAY_HPP
 
 
-#include <cstddef>
-#include <utility>
-#include <vector>
-
-#include "core/utilities/indexer.hpp"
+#include "core/def.hpp"
 #include "particle.hpp"
-#include "core/utilities/point/point.hpp"
-#include "core/utilities/cellmap.hpp"
 #include "core/logger.hpp"
+#include "core/utilities/cellmap.hpp"
 #include "core/utilities/box/box.hpp"
 #include "core/utilities/range/range.hpp"
-#include "core/def.hpp"
+
+
+#include <vector>
+#include <cstddef>
+#include <utility>
+
 
 namespace PHARE::core
 {
@@ -295,10 +295,9 @@ namespace core
         {
         }
 
-        template<typename Container_int, typename Container_double>
-        ContiguousParticles(Container_int&& _iCell, Container_double&& _delta,
-                            Container_double&& _weight, Container_double&& _charge,
-                            Container_double&& _v)
+        template<typename ICell, typename Delta, typename Weight, typename Charge, typename V>
+        ContiguousParticles(ICell&& _iCell, Delta&& _delta, Weight&& _weight, Charge&& _charge,
+                            V&& _v)
             : iCell{_iCell}
             , delta{_delta}
             , weight{_weight}
@@ -319,10 +318,10 @@ namespace core
         NO_DISCARD Return _to(std::size_t i)
         {
             return {
-                *const_cast<double*>(weight.data() + i),     //
-                *const_cast<double*>(charge.data() + i),     //
-                *_array_cast<dim>(iCell.data() + (dim * i)), //
-                *_array_cast<dim>(delta.data() + (dim * i)), //
+                *const_cast<floater_t<3>*>(weight.data() + i), //
+                *const_cast<floater_t<2>*>(charge.data() + i), //
+                *_array_cast<dim>(iCell.data() + (dim * i)),   //
+                *_array_cast<dim>(delta.data() + (dim * i)),   //
                 *_array_cast<3>(v.data() + (3 * i)),
             };
         }
@@ -374,8 +373,10 @@ namespace core
         NO_DISCARD auto cend() const { return iterator(this); }
 
         container_t<int> iCell;
-        container_t<double> delta;
-        container_t<double> weight, charge, v;
+        container_t<floater_t<0>> delta;
+        container_t<floater_t<3>> weight;
+        container_t<floater_t<2>> charge;
+        container_t<floater_t<1>> v;
     };
 
 

src/core/def.hpp

@@ -1,6 +1,8 @@
 #ifndef PHARE_CORE_DEF_HPP
 #define PHARE_CORE_DEF_HPP
 
+#include <bitset>
+#include <cstdint>
 #include <type_traits>
 
 #define NO_DISCARD [[nodiscard]]
@@ -44,6 +46,35 @@ NO_DISCARD bool isSettable(auto const&... args)
     return (check(args) && ...);
 }
 
+
+#ifndef PHARE_DOUBLE_BITSET
+constexpr std::bitset<5> doubles{0b11110}; // index 0 starts on the right in binary
+#else                                      // PHARE_DOUBLE_BITSET
+constexpr std::bitset<5> doubles{PHARE_DOUBLE_BITSET};
+#endif
+
+template<std::uint8_t i>
+bool constexpr is_double()
+{
+    // 0 = particle delta
+    // 1 = particle v
+    // 2 = particle charge
+    // 3 = particle weight
+    // 4 = fields
+
+    return doubles[i] == 1;
+}
+
+template<std::uint8_t i>
+struct Floater
+{
+    using value_type = std::conditional_t<is_double<i>(), double, float>;
+};
+
+template<std::uint8_t i>
+using floater_t = Floater<i>::value_type;
+
+
 } // namespace PHARE::core
 
 #endif // PHARE_CORE_DEF_HPP