Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

try float32 #943

Draft
wants to merge 7 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions res/cmake/dep/highfive.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -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)
PhilipDeegan marked this conversation as resolved.
Show resolved Hide resolved

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
Expand Down
68 changes: 40 additions & 28 deletions src/amr/data/field/field_data.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,10 @@
#define PHARE_SRC_AMR_FIELD_FIELD_DATA_HPP



#include "core/def.hpp"
#include "core/def/phare_mpi.hpp"


#include <SAMRAI/hier/PatchData.h>
#include <SAMRAI/tbox/MemoryUtilities.h>
#include <utility>
Expand Down Expand Up @@ -87,8 +88,16 @@ namespace amr
Super::getFromRestart(restart_db);

assert(field.vector().size() > 0);
restart_db->getDoubleArray("field_" + field.name(), field.vector().data(),
field.vector().size()); // do not reallocate!
if constexpr (std::is_same_v<core::floater_t<4>, double>)
{
restart_db->getDoubleArray("field_" + field.name(), field.vector().data(),
field.vector().size()); // do not reallocate!
}
else
{
restart_db->getFloatArray("field_" + field.name(), field.vector().data(),
field.vector().size()); // do not reallocate!
}
}

void putToRestart(std::shared_ptr<SAMRAI::tbox::Database> const& restart_db) const override
Expand Down Expand Up @@ -471,12 +480,12 @@ namespace amr
void copyImpl(SAMRAI::hier::Box const& localSourceBox, Grid_t const& source,
SAMRAI::hier::Box const& localDestinationBox, Grid_t& destination) const
{
std::uint32_t xSourceStart = static_cast<std::uint32_t>(localSourceBox.lower(0));
std::uint32_t xDestinationStart
std::uint32_t const xSourceStart = static_cast<std::uint32_t>(localSourceBox.lower(0));
std::uint32_t const xDestinationStart
= static_cast<std::uint32_t>(localDestinationBox.lower(0));

std::uint32_t xSourceEnd = static_cast<std::uint32_t>(localSourceBox.upper(0));
std::uint32_t xDestinationEnd
std::uint32_t const xSourceEnd = static_cast<std::uint32_t>(localSourceBox.upper(0));
std::uint32_t const xDestinationEnd
= static_cast<std::uint32_t>(localDestinationBox.upper(0));

for (std::uint32_t xSource = xSourceStart, xDestination = xDestinationStart;
Expand All @@ -489,11 +498,12 @@ namespace amr



void packImpl(std::vector<double>& buffer, Grid_t const& source,
template<typename T>
void packImpl(std::vector<T>& buffer, Grid_t const& source,
SAMRAI::hier::Box const& overlap, SAMRAI::hier::Box const& sourceBox) const
{
int xStart = overlap.lower(0) - sourceBox.lower(0);
int xEnd = overlap.upper(0) - sourceBox.lower(0);
int const xStart = overlap.lower(0) - sourceBox.lower(0);
int const xEnd = overlap.upper(0) - sourceBox.lower(0);

for (int xi = xStart; xi <= xEnd; ++xi)
{
Expand All @@ -502,13 +512,13 @@ namespace amr
}



void unpackImpl(std::size_t& seek, std::vector<double> const& buffer, Grid_t& source,
template<typename T>
void unpackImpl(std::size_t& seek, std::vector<T> const& buffer, Grid_t& source,
SAMRAI::hier::Box const& overlap,
SAMRAI::hier::Box const& destination) const
{
int xStart = overlap.lower(0) - destination.lower(0);
int xEnd = overlap.upper(0) - destination.lower(0);
int const xStart = overlap.lower(0) - destination.lower(0);
int const xEnd = overlap.upper(0) - destination.lower(0);

for (int xi = xStart; xi <= xEnd; ++xi)
{
Expand Down Expand Up @@ -559,16 +569,16 @@ namespace amr




void packImpl(std::vector<double>& buffer, Grid_t const& source,
template<typename T>
void packImpl(std::vector<T>& buffer, Grid_t const& source,
SAMRAI::hier::Box const& overlap, SAMRAI::hier::Box const& destination) const

{
int xStart = overlap.lower(0) - destination.lower(0);
int xEnd = overlap.upper(0) - destination.lower(0);
int const xStart = overlap.lower(0) - destination.lower(0);
int const xEnd = overlap.upper(0) - destination.lower(0);

int yStart = overlap.lower(1) - destination.lower(1);
int yEnd = overlap.upper(1) - destination.lower(1);
int const yStart = overlap.lower(1) - destination.lower(1);
int const yEnd = overlap.upper(1) - destination.lower(1);

for (int xi = xStart; xi <= xEnd; ++xi)
{
Expand All @@ -582,15 +592,16 @@ namespace amr



void unpackImpl(std::size_t& seek, std::vector<double> const& buffer, Grid_t& source,
template<typename T>
void unpackImpl(std::size_t& seek, std::vector<T> const& buffer, Grid_t& source,
SAMRAI::hier::Box const& overlap,
SAMRAI::hier::Box const& destination) const
{
int xStart = overlap.lower(0) - destination.lower(0);
int xEnd = overlap.upper(0) - destination.lower(0);
int const xStart = overlap.lower(0) - destination.lower(0);
int const xEnd = overlap.upper(0) - destination.lower(0);

int yStart = overlap.lower(1) - destination.lower(1);
int yEnd = overlap.upper(1) - destination.lower(1);
int const yStart = overlap.lower(1) - destination.lower(1);
int const yEnd = overlap.upper(1) - destination.lower(1);

for (int xi = xStart; xi <= xEnd; ++xi)
{
Expand Down Expand Up @@ -658,8 +669,8 @@ namespace amr




void packImpl(std::vector<double>& buffer, Grid_t const& source,
template<typename T>
void packImpl(std::vector<T>& buffer, Grid_t const& source,
SAMRAI::hier::Box const& overlap, SAMRAI::hier::Box const& destination) const
{
int xStart = overlap.lower(0) - destination.lower(0);
Expand All @@ -686,7 +697,8 @@ namespace amr



void unpackImpl(std::size_t& seek, std::vector<double> const& buffer, Grid_t& source,
template<typename T>
void unpackImpl(std::size_t& seek, std::vector<T> const& buffer, Grid_t& source,
SAMRAI::hier::Box const& overlap,
SAMRAI::hier::Box const& destination) const
{
Expand Down
1 change: 0 additions & 1 deletion src/core/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
)
Expand Down
32 changes: 13 additions & 19 deletions src/core/data/grid/gridlayout.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,9 @@ namespace core
* @param nbrCells is the number of physical cells of the grid
* @param origin is the point of origin in physical units of the origin of the grid
*/
GridLayout(std::array<double, dimension> const& meshSize,
GridLayout(std::array<floater_t<4>, dimension> const& meshSize,
std::array<std::uint32_t, dimension> const& nbrCells,
Point<double, dimension> const& origin,
Point<floater_t<4>, dimension> const& origin,
Box<int, dimension> AMRBox = Box<int, dimension>{}, int level_number = 0)
: meshSize_{meshSize}
, origin_{origin}
Expand Down Expand Up @@ -155,31 +155,25 @@ namespace core
* @brief origin return the lower point of the grid described by the GridLayout
* in physical coordinates
*/
NO_DISCARD Point<double, dimension> origin() const noexcept { return origin_; }
NO_DISCARD auto origin() const noexcept { return origin_; }



/**
* @brief returns the mesh size in the 'dim' dimensions
*/
NO_DISCARD std::array<double, dimension> const& meshSize() const noexcept
{
return meshSize_;
}
NO_DISCARD auto const& meshSize() const noexcept { return meshSize_; }



NO_DISCARD double inverseMeshSize(Direction direction) const noexcept
NO_DISCARD auto inverseMeshSize(Direction direction) const noexcept
{
return inverseMeshSize_[static_cast<std::uint32_t>(direction)];
}



NO_DISCARD std::array<double, dimension> inverseMeshSize() const noexcept
{
return inverseMeshSize_;
}
NO_DISCARD auto inverseMeshSize() const noexcept { return inverseMeshSize_; }



Expand Down Expand Up @@ -520,28 +514,28 @@ namespace core
* of a multidimensional index that is cell-centered.
*/
template<typename... Indexes>
NO_DISCARD Point<double, dimension> cellCenteredCoordinates(Indexes... index) const
NO_DISCARD auto cellCenteredCoordinates(Indexes... index) const
{
static_assert(sizeof...(Indexes) == dimension,
"Error dimension does not match number of arguments");

std::uint32_t constexpr iPrimal = static_cast<std::uint32_t>(QtyCentering::primal);

constexpr double halfCell = 0.5;
constexpr floater_t<4> halfCell = 0.5;
// A shift of +dx/2, +dy/2, +dz/2 is necessary to get the
// cell center physical coordinates,
// because this point is located on the dual mesh

Point<std::uint32_t, dimension> coord(index...);

Point<double, dimension> physicalPosition;
Point<floater_t<4>, dimension> physicalPosition;

for (std::size_t iDir = 0; iDir < dimension; ++iDir)
{
auto iStart = physicalStartIndexTable_[iPrimal][iDir];

physicalPosition[iDir]
= (static_cast<double>(coord[iDir] - iStart) + halfCell) * meshSize_[iDir]
= (static_cast<floater_t<4>>(coord[iDir] - iStart) + halfCell) * meshSize_[iDir]
+ origin_[iDir];
}

Expand Down Expand Up @@ -1497,10 +1491,10 @@ namespace core



std::array<double, dimension> meshSize_;
Point<double, dimension> origin_;
std::array<floater_t<4>, dimension> meshSize_;
Point<floater_t<4>, dimension> origin_;
std::array<std::uint32_t, dimension> nbrPhysicalCells_;
std::array<double, dimension> inverseMeshSize_;
std::array<floater_t<4>, dimension> inverseMeshSize_;
static constexpr gridDataT data{};

// stores key indices in each direction (3) for primal and dual nodes (2)
Expand Down
4 changes: 2 additions & 2 deletions src/core/data/grid/gridlayoutdefs.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,14 +20,14 @@ namespace core
template<std::size_t dim>
struct WeightPoint
{
constexpr WeightPoint(Point<int, dim> point, double _coef)
constexpr WeightPoint(Point<int, dim> const point, floater_t<4> const _coef)
: indexes{std::move(point)}
, coef{_coef}
{
}

Point<int, dim> indexes;
double coef;
floater_t<4> coef;
};


Expand Down
Original file line number Diff line number Diff line change
@@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
{
Expand All @@ -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);
Expand Down
Loading
Loading