diff --git a/res/cmake/def.cmake b/res/cmake/def.cmake index 684425235..e2a024bbe 100644 --- a/res/cmake/def.cmake +++ b/res/cmake/def.cmake @@ -69,7 +69,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 -Wuninitialized ) + set (_Werr ${_Werr} -Wuninitialized ) # -Wdouble-promotion if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") set (_Werr ${_Werr} -Wno-gnu-zero-variadic-macro-arguments) diff --git a/res/cmake/dep/highfive.cmake b/res/cmake/dep/highfive.cmake index b74081c68..7215a6240 100644 --- a/res/cmake/dep/highfive.cmake +++ b/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 diff --git a/src/amr/data/field/coarsening/coarsen_weighter.hpp b/src/amr/data/field/coarsening/coarsen_weighter.hpp index 0639ff0c8..2720773d4 100644 --- a/src/amr/data/field/coarsening/coarsen_weighter.hpp +++ b/src/amr/data/field/coarsening/coarsen_weighter.hpp @@ -26,12 +26,12 @@ namespace amr computeWeights_(nbrPoints); } - std::vector const& weights() const { return weights_; } + std::vector> const& weights() const { return weights_; } private: - std::vector weights_; + std::vector> weights_; - double findX_(std::size_t nbrPoints) const; + core::floater_t<4> findX_(std::size_t nbrPoints) const; void computeWeights_(std::size_t nbrPoints); }; diff --git a/src/amr/data/field/coarsening/default_field_coarsener.hpp b/src/amr/data/field/coarsening/default_field_coarsener.hpp index ff1356d7f..aa5f14364 100644 --- a/src/amr/data/field/coarsening/default_field_coarsener.hpp +++ b/src/amr/data/field/coarsening/default_field_coarsener.hpp @@ -11,6 +11,7 @@ #include "amr/data/field/coarsening/field_coarsen_index_weight.hpp" #include "amr/resources_manager/amr_utils.hpp" +#include "core/utilities/types.hpp" #include @@ -67,7 +68,7 @@ namespace amr coarseIndex = AMRToLocal(coarseIndex, destinationBox_); - double coarseValue = 0.; + core::floater_t<4> coarseValue = 0.; @@ -99,7 +100,7 @@ namespace amr for (std::size_t iShiftX = 0; iShiftX < xWeights.size(); ++iShiftX) { - double Yinterp = 0.; + core::floater_t<4> Yinterp = 0.; for (std::size_t iShiftY = 0; iShiftY < yWeights.size(); ++iShiftY) { Yinterp += fineField(xStartIndex + iShiftX, yStartIndex + iShiftY) @@ -125,11 +126,11 @@ namespace amr for (std::size_t iShiftX = 0; iShiftX < xWeights.size(); ++iShiftX) { - double Yinterp = 0.; + core::floater_t<4> Yinterp = 0.; for (std::size_t iShiftY = 0; iShiftY < yWeights.size(); ++iShiftY) { - double Zinterp = 0.; + core::floater_t<4> Zinterp = 0.; for (std::size_t iShiftZ = 0; iShiftZ < zWeights.size(); ++iShiftZ) { Zinterp += fineField(xStartIndex + iShiftX, yStartIndex + iShiftY, diff --git a/src/amr/data/field/coarsening/field_coarsen.cpp b/src/amr/data/field/coarsening/field_coarsen.cpp index f5b111afc..f796c5bc6 100644 --- a/src/amr/data/field/coarsening/field_coarsen.cpp +++ b/src/amr/data/field/coarsening/field_coarsen.cpp @@ -2,23 +2,23 @@ namespace PHARE::amr { -NO_DISCARD double CoarsenWeighter::findX_(std::size_t nbrPoints) const +NO_DISCARD core::floater_t<4> CoarsenWeighter::findX_(std::size_t nbrPoints) const { - double x = 0.; + core::floater_t<4> x = 0.; if (nbrPoints % 2 != 0) { - x = 1.; + x = 1.f; for (std::size_t i = 1; i <= (nbrPoints - 1) / 2; ++i) { - x += 2 * 1. / (i + 1); + x += 2 * 1.f / (i + 1); } } else { for (std::size_t i = 1; i <= nbrPoints / 2; ++i) { - x += 2 * 1. / i; + x += 2 * 1.f / i; } } @@ -32,7 +32,7 @@ void CoarsenWeighter::computeWeights_(std::size_t nbrPoints) { weights_.resize(nbrPoints); - auto x = findX_(nbrPoints); + auto const x = findX_(nbrPoints); if (nbrPoints % 2 != 0) @@ -42,13 +42,13 @@ void CoarsenWeighter::computeWeights_(std::size_t nbrPoints) auto const halfNumberOfPointsLeft = nbrPoints / 2; // half of the points needed besides the one on the middle - weights_[halfIndex] = 1. / x; + weights_[halfIndex] = 1.f / x; for (std::size_t i = 1; i <= halfNumberOfPointsLeft; ++i) { - double factor = static_cast(i + 1); - weights_[halfIndex - i] = 1. / (factor * x); - weights_[halfIndex + i] = 1. / (factor * x); + core::floater_t<4> factor = static_cast>(i + 1); + weights_[halfIndex - i] = 1.f / (factor * x); + weights_[halfIndex + i] = 1.f / (factor * x); } } @@ -60,17 +60,17 @@ void CoarsenWeighter::computeWeights_(std::size_t nbrPoints) auto const halfIndexRight = nbrPoints / 2; auto const halfIndexLeft = halfIndexRight - 1; - weights_[halfIndexRight] = 1. / x; - weights_[halfIndexLeft] = 1. / x; + weights_[halfIndexRight] = 1.f / x; + weights_[halfIndexLeft] = 1.f / x; auto const halfNumberOfPointsLeft = (nbrPoints / 2) - 1; for (std::size_t i = 1; i <= halfNumberOfPointsLeft; ++i) { - double factor = static_cast(i + 1); + core::floater_t<4> factor = static_cast>(i + 1); - weights_[halfIndexRight + i] = 1. / (factor * x); - weights_[halfIndexLeft - i] = 1. / (factor * x); + weights_[halfIndexRight + i] = 1.f / (factor * x); + weights_[halfIndexLeft - i] = 1.f / (factor * x); } } } diff --git a/src/amr/data/field/coarsening/field_coarsen_index_weight.hpp b/src/amr/data/field/coarsening/field_coarsen_index_weight.hpp index 79a3b88f8..9d4505acf 100644 --- a/src/amr/data/field/coarsening/field_coarsen_index_weight.hpp +++ b/src/amr/data/field/coarsening/field_coarsen_index_weight.hpp @@ -99,7 +99,7 @@ namespace amr } - NO_DISCARD std::vector const& weights(core::Direction dir) const + NO_DISCARD std::vector> const& weights(core::Direction dir) const { return weighters_[static_cast(dir)].weights(); } diff --git a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp index 39d816413..33ca679d3 100644 --- a/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp +++ b/src/amr/data/field/coarsening/magnetic_field_coarsener.hpp @@ -87,7 +87,8 @@ class MagneticFieldCoarsener else if (centering_[dirX] == core::QtyCentering::dual) // by and bz { coarseField(coarseIndex[dirX]) - = 0.5 * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX])); + = 0.5f + * (fineField(fineStartIndex[dirX] + 1) + fineField(fineStartIndex[dirX])); } } @@ -97,7 +98,7 @@ class MagneticFieldCoarsener and centering_[dirY] == core::QtyCentering::dual) { coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.5 + = 0.5f * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1)); } @@ -105,7 +106,7 @@ class MagneticFieldCoarsener and centering_[dirY] == core::QtyCentering::primal) { coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.5 + = 0.5f * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY])); } @@ -113,7 +114,7 @@ class MagneticFieldCoarsener and centering_[dirY] == core::QtyCentering::dual) { coarseField(coarseIndex[dirX], coarseIndex[dirY]) - = 0.25 + = 0.25f * (fineField(fineStartIndex[dirX], fineStartIndex[dirY]) + fineField(fineStartIndex[dirX] + 1, fineStartIndex[dirY]) + fineField(fineStartIndex[dirX], fineStartIndex[dirY] + 1) diff --git a/src/amr/data/field/field_data.hpp b/src/amr/data/field/field_data.hpp index 385393be7..3f25e5291 100644 --- a/src/amr/data/field/field_data.hpp +++ b/src/amr/data/field/field_data.hpp @@ -2,9 +2,10 @@ #define PHARE_SRC_AMR_FIELD_FIELD_DATA_HPP - +#include "core/def.hpp" #include "core/def/phare_mpi.hpp" + #include #include #include @@ -40,6 +41,7 @@ namespace amr typename PhysicalQuantity = decltype(std::declval().physicalQuantity())> class FieldData : public SAMRAI::hier::PatchData { + static_assert(std::is_same_v>); using Super = SAMRAI::hier::PatchData; public: @@ -62,18 +64,19 @@ namespace amr { } // - [[deprecated]] FieldData(SAMRAI::hier::Box const& domain, - SAMRAI::hier::IntVector const& ghost, std::string name, - std::array const& dl, - std::array const& nbrCells, - core::Point const& origin, PhysicalQuantity qty) + // [[deprecated]] FieldData(SAMRAI::hier::Box const& domain, + // SAMRAI::hier::IntVector const& ghost, std::string name, + // std::array const& dl, + // std::array const& nbrCells, + // core::Point const& origin, PhysicalQuantity + // qty) - : SAMRAI::hier::PatchData(domain, ghost) - , gridLayout{dl, nbrCells, origin} - , field(name, qty, gridLayout.allocSize(qty)) - , quantity_{qty} - { - } + // : SAMRAI::hier::PatchData(domain, ghost) + // , gridLayout{dl, nbrCells, origin} + // , field(name, qty, gridLayout.allocSize(qty)) + // , quantity_{qty} + // { + // } FieldData() = delete; FieldData(FieldData const&) = delete; @@ -87,8 +90,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, 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 const& restart_db) const override @@ -471,12 +482,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(localSourceBox.lower(0)); - std::uint32_t xDestinationStart + std::uint32_t const xSourceStart = static_cast(localSourceBox.lower(0)); + std::uint32_t const xDestinationStart = static_cast(localDestinationBox.lower(0)); - std::uint32_t xSourceEnd = static_cast(localSourceBox.upper(0)); - std::uint32_t xDestinationEnd + std::uint32_t const xSourceEnd = static_cast(localSourceBox.upper(0)); + std::uint32_t const xDestinationEnd = static_cast(localDestinationBox.upper(0)); for (std::uint32_t xSource = xSourceStart, xDestination = xDestinationStart; @@ -489,11 +500,12 @@ namespace amr - void packImpl(std::vector& buffer, Grid_t const& source, + template + void packImpl(std::vector& 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) { @@ -502,13 +514,13 @@ namespace amr } - - void unpackImpl(std::size_t& seek, std::vector const& buffer, Grid_t& source, + template + void unpackImpl(std::size_t& seek, std::vector 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) { @@ -559,16 +571,16 @@ namespace amr - - void packImpl(std::vector& buffer, Grid_t const& source, + template + void packImpl(std::vector& 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) { @@ -582,15 +594,16 @@ namespace amr - void unpackImpl(std::size_t& seek, std::vector const& buffer, Grid_t& source, + template + void unpackImpl(std::size_t& seek, std::vector 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) { @@ -658,8 +671,8 @@ namespace amr - - void packImpl(std::vector& buffer, Grid_t const& source, + template + void packImpl(std::vector& buffer, Grid_t const& source, SAMRAI::hier::Box const& overlap, SAMRAI::hier::Box const& destination) const { int xStart = overlap.lower(0) - destination.lower(0); @@ -686,7 +699,8 @@ namespace amr - void unpackImpl(std::size_t& seek, std::vector const& buffer, Grid_t& source, + template + void unpackImpl(std::size_t& seek, std::vector const& buffer, Grid_t& source, SAMRAI::hier::Box const& overlap, SAMRAI::hier::Box const& destination) const { diff --git a/src/amr/data/field/field_data_factory.hpp b/src/amr/data/field/field_data_factory.hpp index b7c7e0238..202d2da16 100644 --- a/src/amr/data/field/field_data_factory.hpp +++ b/src/amr/data/field/field_data_factory.hpp @@ -35,7 +35,7 @@ namespace amr FieldDataFactory(bool fineBoundaryRepresentsVariable, bool dataLivesOnPatchBorder, std::string const& name, PhysicalQuantity qty) : SAMRAI::hier::PatchDataFactory( - SAMRAI::hier::IntVector{SAMRAI::tbox::Dimension(dimension), n_ghosts}) + SAMRAI::hier::IntVector{SAMRAI::tbox::Dimension(dimension), n_ghosts}) , fineBoundaryRepresentsVariable_{fineBoundaryRepresentsVariable} , dataLivesOnPatchBorder_{dataLivesOnPatchBorder} , quantity_{qty} @@ -88,9 +88,9 @@ namespace amr // for the gridlayout, and also we give the box to the FieldGeometry, so that // it can use it to get the final box representation. - std::array dl; + std::array, dimension> dl; std::array nbCell; - core::Point origin; + core::Point, dimension> origin; for (std::size_t iDim = 0; iDim < dimension; ++iDim) { @@ -116,9 +116,9 @@ namespace amr // TODO: this calculus assumes that we don't need more memory than // alignedMemory(nx*ny*nz*sizeof(double)) + alignedMemory(baseSize) - std::array dl; + std::array, dimension> dl; std::array nbCell; - core::Point origin; + core::Point, dimension> origin; for (std::size_t iDim = 0; iDim < dimension; ++iDim) { diff --git a/src/amr/data/field/refine/electric_field_refiner.hpp b/src/amr/data/field/refine/electric_field_refiner.hpp index 57698b350..24dfaeb13 100644 --- a/src/amr/data/field/refine/electric_field_refiner.hpp +++ b/src/amr/data/field/refine/electric_field_refiner.hpp @@ -129,7 +129,7 @@ class ElectricFieldRefiner // we're on a fine edge in between two coarse edges // we take the average fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); + = 0.5f * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); } } // Ey @@ -150,7 +150,7 @@ class ElectricFieldRefiner // we're on a fine edge in between two coarse ones // we take the average fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); + = 0.5f * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); } } // and this is now Ez @@ -163,13 +163,13 @@ class ElectricFieldRefiner } else if (onCoarseXFace_(fineIndex)) fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); + = 0.5f * (coarseField(ilcx, ilcy) + coarseField(ilcx, ilcy + 1)); else if (onCoarseYFace_(fineIndex)) fineField(ilfx, ilfy) - = 0.5 * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); + = 0.5f * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy)); else fineField(ilfx, ilfy) - = 0.25 + = 0.25f * (coarseField(ilcx, ilcy) + coarseField(ilcx + 1, ilcy) + coarseField(ilcx, ilcy + 1) + coarseField(ilcx + 1, ilcy + 1)); } @@ -208,7 +208,7 @@ class ElectricFieldRefiner else if (onCoarseYFace_(fineIndex)) { fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); } // we share a Z face but not the Y face // we must be one of the 2 X fine edges on a Z face @@ -216,15 +216,15 @@ class ElectricFieldRefiner else if (onCoarseZFace_(fineIndex)) { fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); } else { // we don't share any face thus we're on one of the 2 middle X edges // we take the average of the 4 surrounding X averages fineField(ilfx, ilfy, ilfz) - = 0.25 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)) - + 0.25 + = 0.25f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)) + + 0.25f * (coarseField(ilcx, ilcy, ilcz + 1) + coarseField(ilcx, ilcy + 1, ilcz + 1)); } @@ -251,7 +251,7 @@ class ElectricFieldRefiner // at z and z+dz // take the average of these 2 coarse value fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1)); } // we're on a Z coarse face, but not on a X coarse face // we thus must be one of the 2 Y edges on a Z face @@ -259,7 +259,7 @@ class ElectricFieldRefiner else if (onCoarseZFace_(fineIndex)) { fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); } // now we're not on any of the coarse faces // so we must be one of the two Y edge in the middle of the cell @@ -267,7 +267,7 @@ class ElectricFieldRefiner else { fineField(ilfx, ilfy, ilfz) - = 0.25 + = 0.25f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz) + coarseField(ilcx, ilcy, ilcz + 1) + coarseField(ilcx + 1, ilcy, ilcz + 1)); @@ -290,7 +290,7 @@ class ElectricFieldRefiner else if (onCoarseXFace_(fineIndex)) { fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz)); } // here we're on a coarse Y face, but not a X face // we must be 1 of the 2 Z edges on a Y face @@ -298,7 +298,7 @@ class ElectricFieldRefiner else if (onCoarseYFace_(fineIndex)) { fineField(ilfx, ilfy, ilfz) - = 0.5 * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); + = 0.5f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz)); } // we're not on any coarse face thus must be one of the 2 Z edges // in the middle of the coarse cell @@ -306,7 +306,7 @@ class ElectricFieldRefiner else { fineField(ilfx, ilfy, ilfz) - = 0.25 + = 0.25f * (coarseField(ilcx, ilcy, ilcz) + coarseField(ilcx + 1, ilcy, ilcz) + coarseField(ilcx, ilcy + 1, ilcz + 1) + coarseField(ilcx + 1, ilcy + 1, ilcz)); diff --git a/src/amr/data/field/refine/field_linear_refine.hpp b/src/amr/data/field/refine/field_linear_refine.hpp index 5cb820f45..8d8f0d3a7 100644 --- a/src/amr/data/field/refine/field_linear_refine.hpp +++ b/src/amr/data/field/refine/field_linear_refine.hpp @@ -70,22 +70,22 @@ namespace amr auto coarseIndex{fineIndex}; // here we perform the floating point division, and then we truncate to integer - coarseIndex[dirX] - = std::floor(static_cast(fineIndex[dirX] + shifts_[dirX]) / ratio_(dirX) - - shifts_[dirX]); + coarseIndex[dirX] = std::floor( + static_cast>(fineIndex[dirX] + shifts_[dirX]) / ratio_(dirX) + - shifts_[dirX]); if constexpr (dimension > 1) { - coarseIndex[dirY] - = std::floor(static_cast(fineIndex[dirY] + shifts_[dirY]) / ratio_(dirY) - - shifts_[dirY]); + coarseIndex[dirY] = std::floor( + static_cast>(fineIndex[dirY] + shifts_[dirY]) / ratio_(dirY) + - shifts_[dirY]); } if constexpr (dimension > 2) { - coarseIndex[dirZ] - = std::floor(static_cast(fineIndex[dirZ] + shifts_[dirZ]) / ratio_(dirZ) - - shifts_[dirZ]); + coarseIndex[dirZ] = std::floor( + static_cast>(fineIndex[dirZ] + shifts_[dirZ]) / ratio_(dirZ) + - shifts_[dirZ]); } return coarseIndex; @@ -127,7 +127,7 @@ namespace amr private: SAMRAI::hier::IntVector const ratio_; std::array weighters_; - core::Point shifts_; + core::Point, dimension> shifts_; }; } // namespace amr diff --git a/src/amr/data/field/refine/field_refiner.hpp b/src/amr/data/field/refine/field_refiner.hpp index 89661c08f..b9bfb4a70 100644 --- a/src/amr/data/field/refine/field_refiner.hpp +++ b/src/amr/data/field/refine/field_refiner.hpp @@ -75,7 +75,7 @@ namespace amr coarseStartIndex = AMRToLocal(coarseStartIndex, coarseBox_); fineIndex = AMRToLocal(fineIndex, fineBox_); - double fieldValue = 0.; + core::floater_t<4> fieldValue = 0.; @@ -110,7 +110,7 @@ namespace amr for (std::size_t iShiftX = 0; iShiftX < xLeftRightWeights.size(); ++iShiftX) { - double Yinterp = 0.; + core::floater_t<4> Yinterp = 0.; for (std::size_t iShiftY = 0; iShiftY < yLeftRightWeights.size(); ++iShiftY) { Yinterp += sourceField(xStartIndex + iShiftX, yStartIndex + iShiftY) @@ -142,10 +142,10 @@ namespace amr for (std::size_t iShiftX = 0; iShiftX < xLeftRightWeights.size(); ++iShiftX) { - double Yinterp = 0.; + core::floater_t<4> Yinterp = 0.; for (std::size_t iShiftY = 0; iShiftY < yLeftRightWeights.size(); ++iShiftY) { - double Zinterp = 0.; + core::floater_t<4> Zinterp = 0.; for (std::size_t iShiftZ = 0; iShiftZ < zLeftRightWeights.size(); ++iShiftZ) { Zinterp += sourceField(xStartIndex + iShiftX, yStartIndex + iShiftY, diff --git a/src/amr/data/field/refine/linear_weighter.cpp b/src/amr/data/field/refine/linear_weighter.cpp index c2aa57d63..6dbf96684 100644 --- a/src/amr/data/field/refine/linear_weighter.cpp +++ b/src/amr/data/field/refine/linear_weighter.cpp @@ -11,7 +11,7 @@ LinearWeighter::LinearWeighter(core::QtyCentering centering, std::size_t ratio) assert(nbrPoints > 1); distFromLeftNode_.resize(nbrPoints); bool isEvenRatio = ratio % 2 == 0; - auto smallCellSize = 1. / ratio; + auto smallCellSize = 1.f / ratio; std::iota(std::begin(distFromLeftNode_), std::end(distFromLeftNode_), 0); @@ -20,8 +20,9 @@ LinearWeighter::LinearWeighter(core::QtyCentering centering, std::size_t ratio) if (centering == core::QtyCentering::primal) { std::transform(std::begin(distFromLeftNode_), std::end(distFromLeftNode_), - std::begin(distFromLeftNode_), - [ratio](auto const& v) { return static_cast(v) / ratio; }); + std::begin(distFromLeftNode_), [ratio](auto const& v) { + return static_cast>(v) / ratio; + }); } else { @@ -30,7 +31,7 @@ LinearWeighter::LinearWeighter(core::QtyCentering centering, std::size_t ratio) auto middle = std::begin(distFromLeftNode_) + distFromLeftNode_.size() / 2; std::transform(std::begin(distFromLeftNode_), std::end(distFromLeftNode_), std::begin(distFromLeftNode_), [smallCellSize](auto const& v) { - return (0.5 + static_cast(v)) * smallCellSize; + return (0.5f + static_cast>(v)) * smallCellSize; }); std::rotate(std::begin(distFromLeftNode_), middle, std::end(distFromLeftNode_)); } @@ -40,7 +41,7 @@ LinearWeighter::LinearWeighter(core::QtyCentering centering, std::size_t ratio) auto middle = std::begin(distFromLeftNode_) + distFromLeftNode_.size() / 2 + 1; std::transform(std::begin(distFromLeftNode_), std::end(distFromLeftNode_), std::begin(distFromLeftNode_), [smallCellSize](auto const& v) { - return static_cast(v) * smallCellSize; + return static_cast>(v) * smallCellSize; }); std::rotate(std::begin(distFromLeftNode_), middle, std::end(distFromLeftNode_)); @@ -49,9 +50,8 @@ LinearWeighter::LinearWeighter(core::QtyCentering centering, std::size_t ratio) std::transform(std::begin(distFromLeftNode_), std::end(distFromLeftNode_), - std::back_inserter(weights_), [](auto const& d) { - return std::array{{1. - d, d}}; - }); + std::back_inserter(weights_), + [](auto const& d) { return std::array, 2>{{1.f - d, d}}; }); } } // namespace PHARE::amr diff --git a/src/amr/data/field/refine/linear_weighter.hpp b/src/amr/data/field/refine/linear_weighter.hpp index 597ac4792..a4a654f35 100644 --- a/src/amr/data/field/refine/linear_weighter.hpp +++ b/src/amr/data/field/refine/linear_weighter.hpp @@ -27,17 +27,20 @@ namespace amr class LinearWeighter { public: - using FineIndexWeight = std::array; + using FineIndexWeight = std::array, 2>; using FineIndexWeights = std::vector; LinearWeighter(core::QtyCentering centering, std::size_t ratio); - std::vector const& getUniformDistances() const { return distFromLeftNode_; } + std::vector> const& getUniformDistances() const + { + return distFromLeftNode_; + } FineIndexWeights const& weights() const { return weights_; } private: - std::vector distFromLeftNode_; + std::vector> distFromLeftNode_; FineIndexWeights weights_; }; diff --git a/src/amr/data/field/refine/magnetic_field_refiner.hpp b/src/amr/data/field/refine/magnetic_field_refiner.hpp index 0369c857a..4dd8bbaba 100644 --- a/src/amr/data/field/refine/magnetic_field_refiner.hpp +++ b/src/amr/data/field/refine/magnetic_field_refiner.hpp @@ -76,7 +76,7 @@ class MagneticFieldRefiner else { fineField(locFineIdx[dirX]) - = 0.5 + = 0.5f * (coarseField(locCoarseIdx[dirX]) + coarseField(locCoarseIdx[dirX] + 1)); } } @@ -115,7 +115,7 @@ class MagneticFieldRefiner // we're on a fine X face, take the average of the coarse value // between the two surrounding faces fineField(locFineIdx[dirX], locFineIdx[dirY]) - = 0.5 + = 0.5f * (coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]) + coarseField(locCoarseIdx[dirX] + 1, locCoarseIdx[dirY])); } @@ -136,7 +136,7 @@ class MagneticFieldRefiner // we're on a fine Y face, take the average of the coarse value // between the two surrounding faces fineField(locFineIdx[dirX], locFineIdx[dirY]) - = 0.5 + = 0.5f * (coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY]) + coarseField(locCoarseIdx[dirX], locCoarseIdx[dirY] + 1)); } @@ -176,7 +176,7 @@ class MagneticFieldRefiner // we're on a fine X face, take the average of the coarse value // between the two surrounding faces fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = 0.5 * (coarseField(ix, iy, iz) + coarseField(ix + 1, iy, iz)); + = 0.5f * (coarseField(ix, iy, iz) + coarseField(ix + 1, iy, iz)); } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -196,7 +196,7 @@ class MagneticFieldRefiner // we're on a fine Y face, take the average of the coarse value // between the two surrounding faces fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = 0.5 * (coarseField(ix, iy, iz) + coarseField(ix, iy + 1, iz)); + = 0.5f * (coarseField(ix, iy, iz) + coarseField(ix, iy + 1, iz)); } } else if (centerings_[dirX] == core::QtyCentering::dual @@ -216,7 +216,7 @@ class MagneticFieldRefiner // we're on a fine Z face, take the average of the coarse value // between the two surrounding faces fineField(locFineIdx[dirX], locFineIdx[dirY], locFineIdx[dirZ]) - = 0.5 * (coarseField(ix, iy, iz) + coarseField(ix, iy, iz + 1)); + = 0.5f * (coarseField(ix, iy, iz) + coarseField(ix, iy, iz + 1)); } } } diff --git a/src/amr/data/field/time_interpolate/field_linear_time_interpolate.hpp b/src/amr/data/field/time_interpolate/field_linear_time_interpolate.hpp index 3ec28b75d..9d5c4324b 100644 --- a/src/amr/data/field/time_interpolate/field_linear_time_interpolate.hpp +++ b/src/amr/data/field/time_interpolate/field_linear_time_interpolate.hpp @@ -52,10 +52,10 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator auto const& fieldDataSrcOld = dynamic_cast(srcDataOld); auto const& fieldDataSrcNew = dynamic_cast(srcDataNew); - double const interpTime = fieldDataDest.getTime(); - double const oldTime = fieldDataSrcOld.getTime(); - double const newTime = fieldDataSrcNew.getTime(); - double const alpha = (interpTime - oldTime) / (newTime - oldTime); + core::floater_t<4> const interpTime = fieldDataDest.getTime(); + core::floater_t<4> const oldTime = fieldDataSrcOld.getTime(); + core::floater_t<4> const newTime = fieldDataSrcNew.getTime(); + core::floater_t<4> const alpha = (interpTime - oldTime) / (newTime - oldTime); auto const& fieldSrcOld = fieldDataSrcOld.field; auto const& fieldSrcNew = fieldDataSrcNew.field; @@ -90,7 +90,7 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator for (auto ix = iDestStartX, ixSrc = iSrcStartX; ix <= iDestEndX; ++ix, ++ixSrc) { - fieldDest(ix) = (1. - alpha) * fieldSrcOld(ixSrc) + alpha * fieldSrcNew(ixSrc); + fieldDest(ix) = (1.f - alpha) * fieldSrcOld(ixSrc) + alpha * fieldSrcNew(ixSrc); } } else if constexpr (dim == 2) @@ -107,7 +107,7 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator { for (auto iy = iDestStartY, iySrc = iSrcStartY; iy <= iDestEndY; ++iy, ++iySrc) { - fieldDest(ix, iy) = (1. - alpha) * fieldSrcOld(ixSrc, iySrc) + fieldDest(ix, iy) = (1.f - alpha) * fieldSrcOld(ixSrc, iySrc) + alpha * fieldSrcNew(ixSrc, iySrc); } } @@ -131,7 +131,7 @@ class FieldLinearTimeInterpolate : public SAMRAI::hier::TimeInterpolateOperator { for (auto iz = iDestStartZ, izSrc = iSrcStartZ; iz <= iDestEndZ; ++iz, ++izSrc) { - fieldDest(ix, iy, iz) = (1. - alpha) * fieldSrcOld(ixSrc, iySrc, izSrc) + fieldDest(ix, iy, iz) = (1.f - alpha) * fieldSrcOld(ixSrc, iySrc, izSrc) + alpha * fieldSrcNew(ixSrc, iySrc, izSrc); } } diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp index 0efd3f795..776d9fa48 100644 --- a/src/amr/resources_manager/amr_utils.hpp +++ b/src/amr/resources_manager/amr_utils.hpp @@ -141,9 +141,9 @@ namespace amr // We get geometry information from the patch, such as meshSize, and physical origin auto patchGeom = std::dynamic_pointer_cast( patch.getPatchGeometry()); - core::Point origin; + core::Point, dimension> origin; - std::array dl; + std::array, dimension> dl; if (patchGeom != nullptr) { diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp index b9a016f7d..35eb70b4f 100644 --- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp +++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp @@ -16,16 +16,17 @@ class DefaultHybridTaggerStrategy : public HybridTaggerStrategy using gridlayout_type = typename HybridModel::gridlayout_type; static auto constexpr dimension = HybridModel::dimension; + static constexpr core::floater_t<4> def_threshold_ = 0.1; public: DefaultHybridTaggerStrategy(initializer::PHAREDict const& dict) - : threshold_{cppdict::get_value(dict, "threshold", 0.1)} + : threshold_{cppdict::get_value(dict, "threshold", def_threshold_)} { } void tag(HybridModel& model, gridlayout_type const& layout, int* tags) const override; private: - double threshold_ = 0.1; + core::floater_t<4> threshold_ = def_threshold_; }; template @@ -78,10 +79,10 @@ void DefaultHybridTaggerStrategy::tag(HybridModel& model, std::size_t oneOrZero = doLastCell ? 1 : 0; for (auto iCell = 0u, ix = start_x; iCell < end_x + oneOrZero; ++ix, ++iCell) { - auto Byavg = 0.2 * (By(ix - 2) + By(ix - 1) + By(ix) + By(ix + 1) + By(ix + 2)); - auto Bzavg = 0.2 * (Bz(ix - 2) + Bz(ix - 1) + Bz(ix) + Bz(ix + 1) + Bz(ix + 2)); - auto Byavgp1 = 0.2 * (By(ix - 1) + By(ix) + By(ix + 1) + By(ix + 2) + By(ix + 3)); - auto Bzavgp1 = 0.2 * (Bz(ix - 1) + Bz(ix) + Bz(ix + 1) + Bz(ix + 2) + Bz(ix + 3)); + auto Byavg = 0.2f * (By(ix - 2) + By(ix - 1) + By(ix) + By(ix + 1) + By(ix + 2)); + auto Bzavg = 0.2f * (Bz(ix - 2) + Bz(ix - 1) + Bz(ix) + Bz(ix + 1) + Bz(ix + 2)); + auto Byavgp1 = 0.2f * (By(ix - 1) + By(ix) + By(ix + 1) + By(ix + 2) + By(ix + 3)); + auto Bzavgp1 = 0.2f * (Bz(ix - 1) + Bz(ix) + Bz(ix + 1) + Bz(ix + 2) + Bz(ix + 3)); auto criter_by = std::abs(Byavgp1 - Byavg) / (1 + std::abs(Byavg)); auto criter_bz = std::abs(Bzavgp1 - Bzavg) / (1 + std::abs(Bzavg)); auto criter_b = std::sqrt(criter_by * criter_by + criter_bz * criter_bz); diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 293de3711..ee0ccbb04 100644 --- a/src/core/CMakeLists.txt +++ b/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 ) diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp index a079022da..6a732ac56 100644 --- a/src/core/data/electrons/electrons.hpp +++ b/src/core/data/electrons/electrons.hpp @@ -155,7 +155,7 @@ class IsothermalElectronPressureClosure IsothermalElectronPressureClosure(PHARE::initializer::PHAREDict const& dict, Ions const& ions) : ions_{ions} - , Te_{dict["Te"].template to()} + , Te_{dict["Te"].template to>()} , Pe_{"Pe", HybridQuantity::Scalar::P} { } @@ -218,7 +218,7 @@ class IsothermalElectronPressureClosure private: Ions ions_; - double const Te_ = 0; + floater_t<4> const Te_ = 0; Field Pe_; }; diff --git a/src/core/data/field/initializers/field_user_initializer.hpp b/src/core/data/field/initializers/field_user_initializer.hpp index d9c74275a..ea8a8ee8e 100644 --- a/src/core/data/field/initializers/field_user_initializer.hpp +++ b/src/core/data/field/initializers/field_user_initializer.hpp @@ -22,9 +22,9 @@ class FieldUserFunctionInitializer return gridLayout.fieldNodeCoordinates(field_, gridLayout.origin(), args...); }); - std::shared_ptr> gridPtr // keep grid data alive + std::shared_ptr>> const gridPtr // keep grid data alive = std::apply([&](auto&... args) { return init(args...); }, coords); - Span& grid = *gridPtr; + auto& grid = *gridPtr; for (std::size_t cell_idx = 0; cell_idx < indices.size(); cell_idx++) std::apply([&](auto&... args) { field(args...) = grid[cell_idx]; }, indices[cell_idx]); diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 74297a864..fd43376af 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -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 const& meshSize, + GridLayout(std::array, dimension> const& meshSize, std::array const& nbrCells, - Point const& origin, + Point, dimension> const& origin, Box AMRBox = Box{}, int level_number = 0) : meshSize_{meshSize} , origin_{origin} @@ -135,13 +135,13 @@ namespace core } - inverseMeshSize_[0] = 1. / meshSize_[0]; + inverseMeshSize_[0] = 1.f / meshSize_[0]; if constexpr (dimension > 1) { - inverseMeshSize_[1] = 1. / meshSize_[1]; + inverseMeshSize_[1] = 1.f / meshSize_[1]; if constexpr (dimension > 2) { - inverseMeshSize_[2] = 1. / meshSize_[2]; + inverseMeshSize_[2] = 1.f / meshSize_[2]; } } } @@ -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 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 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(direction)]; } - NO_DISCARD std::array inverseMeshSize() const noexcept - { - return inverseMeshSize_; - } + NO_DISCARD auto inverseMeshSize() const noexcept { return inverseMeshSize_; } @@ -278,7 +272,7 @@ namespace core std::get<2>(vectors).emplace_back(point[2]); }; - auto xyz = tuple_fixed_type, dimension>{}; + auto xyz = tuple_fixed_type>, dimension>{}; for (auto const& indiceTuple : indices) std::apply( @@ -293,10 +287,10 @@ namespace core } - NO_DISCARD double cellVolume() const + NO_DISCARD floater_t<4> cellVolume() const { return std::accumulate(meshSize().begin(), meshSize().end(), 1.0, - std::multiplies()); + std::multiplies>()); } /** @@ -466,31 +460,30 @@ namespace core * associated with a given Field, in physical coordinates. */ template - NO_DISCARD Point - fieldNodeCoordinates(Field_t const& field, Point const& origin, - Indexes... index) const + NO_DISCARD auto fieldNodeCoordinates(Field_t const& field, + Point, dimension> const& origin, + Indexes... index) const { static_assert(sizeof...(Indexes) == dimension, "Error dimension does not match number of arguments"); - std::uint32_t iQuantity = static_cast(field.physicalQuantity()); + std::uint32_t const iQuantity = static_cast(field.physicalQuantity()); constexpr std::uint32_t iDual = static_cast(QtyCentering::dual); constexpr auto& hybridQtyCentering = GridLayoutImpl::hybridQtyCentering_; - Point coord{static_cast(index)...}; + Point const coord{static_cast(index)...}; - Point position; + Point, dimension> position; for (std::size_t iDir = 0; iDir < dimension; ++iDir) { - double halfCell = 0.0; - auto const centering = static_cast(hybridQtyCentering[iQuantity][iDir]); - std::int32_t const iStart = physicalStartIndexTable_[centering][iDir]; + std::int32_t const iStart = physicalStartIndexTable_[centering][iDir]; + floater_t<4> const halfCell = centering == iDual ? .5f : 0.0f; // A shift of +dx/2, +dy/2, +dz/2 is necessary to get the physical // coordinate on the dual mesh @@ -501,13 +494,8 @@ namespace core // if ix is dual then ixStart is dual // if iy is primal then iyStart is primal ... - if (centering == iDual) - { - halfCell = 0.5; - } - position[iDir] - = (static_cast(coord[iDir] - iStart) + halfCell) * meshSize_[iDir] + = (static_cast>(coord[iDir] - iStart) + halfCell) * meshSize_[iDir] + origin[iDir]; } @@ -520,28 +508,28 @@ namespace core * of a multidimensional index that is cell-centered. */ template - NO_DISCARD Point 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(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 coord(index...); - Point physicalPosition; + Point, dimension> physicalPosition; for (std::size_t iDir = 0; iDir < dimension; ++iDir) { auto iStart = physicalStartIndexTable_[iPrimal][iDir]; physicalPosition[iDir] - = (static_cast(coord[iDir] - iStart) + halfCell) * meshSize_[iDir] + = (static_cast>(coord[iDir] - iStart) + halfCell) * meshSize_[iDir] + origin_[iDir]; } @@ -628,7 +616,7 @@ namespace core * on the dimensionality of the GridLayout. */ template - NO_DISCARD auto deriv(Field const& operand, MeshIndex index) + NO_DISCARD floater_t<4> deriv(Field const& operand, MeshIndex index) { auto fieldCentering = centering(operand.physicalQuantity()); using PHARE::core::dirX; @@ -695,8 +683,10 @@ namespace core * on the dimensionality of the GridLayout. */ template - NO_DISCARD auto laplacian(Field const& operand, MeshIndex index) + NO_DISCARD floater_t<4> laplacian(Field const& operand, MeshIndex index) { + floater_t<4> constexpr static _2 = 2.0; + static_assert(Field::dimension == dimension, "field dimension must be equal to gridlayout dimension"); using PHARE::core::dirX; @@ -710,7 +700,7 @@ namespace core auto nextX = operand(index[0] + 1); return inverseMeshSize_[dirX] * inverseMeshSize_[dirX] - * (nextX - 2.0 * hereX + prevX); + * (nextX - _2 * hereX + prevX); } else if constexpr (Field::dimension == 2) @@ -720,14 +710,14 @@ namespace core auto nextX = operand(index[0] + 1, index[1]); auto lapX = inverseMeshSize_[dirX] * inverseMeshSize_[dirX] - * (nextX - 2.0 * hereX + prevX); + * (nextX - _2 * hereX + prevX); auto prevY = operand(index[0], index[1] - 1); auto hereY = operand(index[0], index[1]); auto nextY = operand(index[0], index[1] + 1); auto lapY = inverseMeshSize_[dirY] * inverseMeshSize_[dirY] - * (nextY - 2.0 * hereY + prevY); + * (nextY - _2 * hereY + prevY); return lapX + lapY; } @@ -738,21 +728,21 @@ namespace core auto nextX = operand(index[0] + 1, index[1], index[2]); auto lapX = inverseMeshSize_[dirX] * inverseMeshSize_[dirX] - * (nextX - 2.0 * hereX + prevX); + * (nextX - _2 * hereX + prevX); auto prevY = operand(index[0], index[1] - 1, index[2]); auto hereY = operand(index[0], index[1], index[2]); auto nextY = operand(index[0], index[1] + 1, index[2]); auto lapY = inverseMeshSize_[dirY] * inverseMeshSize_[dirY] - * (nextY - 2.0 * hereY + prevY); + * (nextY - _2 * hereY + prevY); auto prevZ = operand(index[0], index[1], index[2] - 1); auto hereZ = operand(index[0], index[1], index[2]); auto nextZ = operand(index[0], index[1], index[2] + 1); auto lapZ = inverseMeshSize_[dirZ] * inverseMeshSize_[dirZ] - * (nextZ - 2.0 * hereZ + prevZ); + * (nextZ - _2 * hereZ + prevZ); return lapX + lapY + lapZ; } @@ -1497,10 +1487,10 @@ namespace core - std::array meshSize_; - Point origin_; + std::array, dimension> meshSize_; + Point, dimension> origin_; std::array nbrPhysicalCells_; - std::array inverseMeshSize_; + std::array, dimension> inverseMeshSize_; static constexpr gridDataT data{}; // stores key indices in each direction (3) for primal and dual nodes (2) diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp index 23920396e..6186e2b7a 100644 --- a/src/core/data/grid/gridlayoutdefs.hpp +++ b/src/core/data/grid/gridlayoutdefs.hpp @@ -20,14 +20,14 @@ namespace core template struct WeightPoint { - constexpr WeightPoint(Point point, double _coef) + constexpr WeightPoint(Point const point, floater_t<4> const _coef) : indexes{std::move(point)} , coef{_coef} { } Point indexes; - double coef; + floater_t<4> coef; }; diff --git a/src/core/data/ions/ion_population/ion_population.hpp b/src/core/data/ions/ion_population/ion_population.hpp index aa6554e47..7b4b92bc3 100644 --- a/src/core/data/ions/ion_population/ion_population.hpp +++ b/src/core/data/ions/ion_population/ion_population.hpp @@ -28,7 +28,7 @@ namespace core IonPopulation(initializer::PHAREDict const& initializer) : name_{initializer["name"].template to()} - , mass_{initializer["mass"].template to()} + , mass_{static_cast>(initializer["mass"].template to())} , flux_{name_ + "_flux", HybridQuantity::Vector::V} , momentumTensor_{name_ + "_momentumTensor", HybridQuantity::Tensor::M} , rho_{name_ + "_rho", HybridQuantity::Scalar::rho} @@ -38,7 +38,7 @@ namespace core } - NO_DISCARD double mass() const { return mass_; } + NO_DISCARD floater_t<4> mass() const { return mass_; } NO_DISCARD std::string const& name() const { return name_; } @@ -121,7 +121,7 @@ namespace core private: std::string name_; - double mass_; + floater_t<4> mass_; VecField flux_; TensorField momentumTensor_; field_type rho_; diff --git a/src/core/data/ions/ions.hpp b/src/core/data/ions/ions.hpp index 2701ebdac..d2ea54931 100644 --- a/src/core/data/ions/ions.hpp +++ b/src/core/data/ions/ions.hpp @@ -257,8 +257,9 @@ namespace core private: bool allSameMass_() const { + floater_t<4> constexpr diff = 1e-6; return all(populations_, [this](auto const& pop) { // arbitrary small diff - return float_equals(pop.mass(), populations_.front().mass(), /*abs_tol=*/1e-10); + return float_equals(pop.mass(), populations_.front().mass(), /*abs_tol=*/diff); }); } diff --git a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp index 49098af06..8e79e79aa 100644 --- a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp +++ b/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 -#include -#include -#include -#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 const& V, std::array const& Vth, - std::mt19937_64& generator, std::array& 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 basisTransform(std::array, 3> const& basis, - std::array const& vec); +#include +#include +#include -void localMagneticBasis(std::array B, std::array, 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. @@ -40,11 +32,11 @@ class MaxwellianParticleInitializer : public ParticleInitializer const& bulkVelocity, - std::array const& thermalVelocity, double const particleCharge, + std::array const& thermalVelocity, floater_t<4> const particleCharge, std::uint32_t const& nbrParticlesPerCell, std::optional seed = {}, Basis const basis = Basis::Cartesian, std::array const& magneticField = {nullptr, nullptr, nullptr}, - double densityCutOff = 1e-5) + floater_t<4> densityCutOff = 1e-5) : density_{density} , bulkVelocity_{bulkVelocity} , thermalVelocity_{thermalVelocity} @@ -87,8 +79,8 @@ class MaxwellianParticleInitializer : public ParticleInitializer thermalVelocity_; std::array magneticField_; - double particleCharge_; - double densityCutOff_ = 1e-5; + floater_t<4> particleCharge_; + floater_t<4> densityCutOff_ = 1e-5; std::uint32_t nbrParticlePerCell_; Basis basis_; std::optional rngSeed_; @@ -116,19 +108,19 @@ class MaxwellianInitFunctions _B[i] = magneticField[i](coords...); } - NO_DISCARD std::array B() const { return ptrs(_B); } + NO_DISCARD std::array const*, 3> B() const { return ptrs(_B); } NO_DISCARD auto operator()() const { return std::make_tuple(_n->data(), ptrs(_V), ptrs(_Vth)); } private: - NO_DISCARD std::array - ptrs(std::array>, 3> const& v) const + NO_DISCARD std::array const*, 3> + ptrs(std::array>>, 3> const& v) const { return {v[0]->data(), v[1]->data(), v[2]->data()}; } - std::shared_ptr> const _n, _ppc; - std::array>, 3> _B, _V, _Vth; + std::shared_ptr>> const _n, _ppc; + std::array>>, 3> _B, _V, _Vth; }; @@ -148,7 +140,7 @@ void MaxwellianParticleInitializer::loadParticles( }; - auto deltas = [](auto& pos, auto& gen) -> std::array { + auto deltas = [](auto& pos, auto& gen) -> std::array, dimension> { if constexpr (dimension == 1) return {pos(gen)}; if constexpr (dimension == 2) @@ -176,7 +168,7 @@ void MaxwellianParticleInitializer::loadParticles( auto const [n, V, Vth] = fns(); auto randGen = getRNG(rngSeed_); - ParticleDeltaDistribution deltaDistrib; + ParticleDeltaDistribution> deltaDistrib; for (std::size_t flatCellIdx = 0; flatCellIdx < ndCellIndices.size(); ++flatCellIdx) { @@ -186,20 +178,20 @@ void MaxwellianParticleInitializer::loadParticles( auto const cellWeight = n[flatCellIdx] / nbrParticlePerCell_; auto const AMRCellIndex = layout.localToAMR(point(flatCellIdx, ndCellIndices)); auto const iCell = AMRCellIndex.template toArray(); - std::array particleVelocity; - std::array, 3> basis; + std::array, 3> particleVelocity; + std::array, 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); diff --git a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.cpp b/src/core/data/ions/particle_initializers/maxwellian_particle_initializer_base.hpp similarity index 61% rename from src/core/data/ions/particle_initializers/maxwellian_particle_initializer.cpp rename to src/core/data/ions/particle_initializers/maxwellian_particle_initializer_base.hpp index 4d528609f..ee13bf554 100644 --- a/src/core/data/ions/particle_initializers/maxwellian_particle_initializer.cpp +++ b/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 +#include #include -#include "core/utilities/types.hpp" #include "core/def.hpp" +#include "core/utilities/types.hpp" namespace PHARE { namespace core { - void maxwellianVelocity(std::array const& V, std::array const& Vth, - std::mt19937_64& generator, std::array& partVelocity) + + template + void maxwellianVelocity(std::array& 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 - basisTransform(std::array, 3> const& basis, - std::array const& vec) + template + NO_DISCARD std::array basisTransform(std::array, 3> const& basis, + std::array const& vec) { - std::array newVec; + std::array newVec; for (std::uint32_t comp = 0; comp < 3; comp++) { @@ -40,12 +47,16 @@ namespace core } - - void localMagneticBasis(std::array B, std::array, 3>& basis) + template + void localMagneticBasis(std::array, 3>& basis, Bargs const... bargs) { + floater_t<4> constexpr static _1em8 = 1e-8; + + std::array const B{bargs...}; + auto b2 = norm(B); - if (b2 < 1e-8) + if (b2 < _1em8) { basis[0][0] = 1.0; basis[0][1] = 0.0; @@ -85,3 +96,5 @@ namespace core } // namespace core } // namespace PHARE + +#endif /*PHARE_MAXWELLIAN_PARTICLE_INITIALIZER_BASE_HPP*/ diff --git a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp b/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp index 97cfa6b0c..36f7cebb1 100644 --- a/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp +++ b/src/core/data/ions/particle_initializers/particle_initializer_factory.hpp @@ -45,9 +45,14 @@ namespace core auto& vthz = dict["thermal_velocity_z"].template to(); - auto charge = dict["charge"].template to(); - - auto densityCutOff = cppdict::get_value(dict, "density_cut_off", double{1e-16}); + auto charge = dict["charge"].template to>(); + + auto const densityCutOff = [&]() { + if constexpr (std::is_same_v, double>) + return cppdict::get_value(dict, "density_cut_off", double{1e-16}); + else + return cppdict::get_value(dict, "density_cut_off", float{1e-6}); + }(); auto nbrPartPerCell = dict["nbr_part_per_cell"].template to(); diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index 7d52d1d01..98d6fe48e 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -2,16 +2,16 @@ #define PHARE_CORE_DATA_NDARRAY_NDARRAY_VECTOR_HPP #include "core/def.hpp" -#include + +#include "core/utilities/types.hpp" + +// #include #include -#include -#include #include +#include +#include #include -#include - - -#include "core/utilities/types.hpp" +#include namespace PHARE::core diff --git a/src/core/data/particles/particle.hpp b/src/core/data/particles/particle.hpp index 0457865c2..2297d1f69 100644 --- a/src/core/data/particles/particle.hpp +++ b/src/core/data/particles/particle.hpp @@ -20,12 +20,15 @@ namespace PHARE::core template struct ParticleDeltaDistribution { + T constexpr static one = 1; + T constexpr static zro = 0; + template NO_DISCARD T operator()(Generator& generator) { return dist(generator); } - std::uniform_real_distribution dist{0, 1. - std::numeric_limits::epsilon()}; + std::uniform_real_distribution dist{zro, one - std::numeric_limits::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 cell, - std::array a_delta, std::array a_v) + Particle(floater_t<3> a_weight, floater_t<2> a_charge, std::array cell, + std::array, dim> a_delta, std::array, 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 iCell = ConstArray(); - std::array delta = ConstArray(); - std::array v = ConstArray(); + std::array iCell = ConstArray(); + std::array, dim> delta = ConstArray, dim>(); + std::array, 3> v = ConstArray, 3>(); NO_DISCARD bool operator==(Particle 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& iCell; - std::array& delta; - std::array& v; + std::array, dim>& delta; + std::array, 3>& v; }; @@ -121,9 +124,9 @@ inline constexpr auto is_phare_particle_type template typename ParticleA, template typename ParticleB> -NO_DISCARD typename std::enable_if_t< - is_phare_particle_type> and is_phare_particle_type>, - bool> +NO_DISCARD typename std::enable_if_t> + and is_phare_particle_type>, + bool> operator==(ParticleA const& particleA, ParticleB const& particleB) { return particleA.weight == particleB.weight and // diff --git a/src/core/data/particles/particle_array.hpp b/src/core/data/particles/particle_array.hpp index 93acc3066..96ce983e9 100644 --- a/src/core/data/particles/particle_array.hpp +++ b/src/core/data/particles/particle_array.hpp @@ -2,18 +2,18 @@ #define PHARE_CORE_DATA_PARTICLES_PARTICLE_ARRAY_HPP -#include -#include -#include - -#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 +#include +#include + namespace PHARE::core { @@ -295,10 +295,9 @@ namespace core { } - template - ContiguousParticles(Container_int&& _iCell, Container_double&& _delta, - Container_double&& _weight, Container_double&& _charge, - Container_double&& _v) + template + 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(weight.data() + i), // - *const_cast(charge.data() + i), // - *_array_cast(iCell.data() + (dim * i)), // - *_array_cast(delta.data() + (dim * i)), // + *const_cast*>(weight.data() + i), // + *const_cast*>(charge.data() + i), // + *_array_cast(iCell.data() + (dim * i)), // + *_array_cast(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 iCell; - container_t delta; - container_t weight, charge, v; + container_t> delta; + container_t> weight; + container_t> charge; + container_t> v; }; diff --git a/src/core/data/particles/particle_utilities.hpp b/src/core/data/particles/particle_utilities.hpp index cdbc4e097..177644956 100644 --- a/src/core/data/particles/particle_utilities.hpp +++ b/src/core/data/particles/particle_utilities.hpp @@ -19,7 +19,7 @@ template NO_DISCARD auto positionAsPoint(Particle const& particle, GridLayout const& layout) { - Point position; + Point, GridLayout::dimension> position; auto origin = layout.origin(); auto startIndexes = layout.physicalStartIndex(QtyCentering::primal); auto meshSize = layout.meshSize(); diff --git a/src/core/def.hpp b/src/core/def.hpp index b741cdb43..fbaef32a6 100644 --- a/src/core/def.hpp +++ b/src/core/def.hpp @@ -1,6 +1,8 @@ #ifndef PHARE_CORE_DEF_HPP #define PHARE_CORE_DEF_HPP +#include +#include #include #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{0b00000}; // index 0 starts on the right in binary +#else // PHARE_DOUBLE_BITSET +constexpr std::bitset<5> doubles{PHARE_DOUBLE_BITSET}; +#endif + +template +bool constexpr is_double() +{ + // 0 = particle delta + // 1 = particle v + // 2 = particle charge + // 3 = particle weight + // 4 = fields + + return 0; // doubles[i] == 1; +} + +template +struct Floater +{ + using value_type = std::conditional_t(), double, float>; +}; + +template +using floater_t = Floater::value_type; + + } // namespace PHARE::core #endif // PHARE_CORE_DEF_HPP diff --git a/src/core/def/phare_mpi.hpp b/src/core/def/phare_mpi.hpp index 4588dde90..0c0bdb843 100644 --- a/src/core/def/phare_mpi.hpp +++ b/src/core/def/phare_mpi.hpp @@ -14,7 +14,7 @@ // clang-format off DISABLE_WARNING(cast-function-type, bad-function-cast, 42) -#include "mpi.h" +#include "mpi.h" // IWYU pragma: private, include "core/def/phare_mpi.hpp" ENABLE_WARNING(cast-function-type, bad-function-cast, 42) // clang-format on diff --git a/src/core/hybrid/hybrid_quantities.hpp b/src/core/hybrid/hybrid_quantities.hpp index ab55eff9d..add60a17e 100644 --- a/src/core/hybrid/hybrid_quantities.hpp +++ b/src/core/hybrid/hybrid_quantities.hpp @@ -69,19 +69,6 @@ class HybridQuantity // no condition, for now there's only then momentum tensor M return {{Scalar::Mxx, Scalar::Mxy, Scalar::Mxz, Scalar::Myy, Scalar::Myz, Scalar::Mzz}}; } - - NO_DISCARD static constexpr auto B_items() - { - auto const& [Bx, By, Bz] = B(); - return std::make_tuple(std::make_pair("Bx", Bx), std::make_pair("By", By), - std::make_pair("Bz", Bz)); - } - NO_DISCARD static constexpr auto E_items() - { - auto const& [Ex, Ey, Ez] = E(); - return std::make_tuple(std::make_pair("Ex", Ex), std::make_pair("Ey", Ey), - std::make_pair("Ez", Ez)); - } }; } // namespace PHARE::core diff --git a/src/core/numerics/faraday/faraday.hpp b/src/core/numerics/faraday/faraday.hpp index 1b9b4e2b5..de283776f 100644 --- a/src/core/numerics/faraday/faraday.hpp +++ b/src/core/numerics/faraday/faraday.hpp @@ -18,7 +18,7 @@ class Faraday : public LayoutHolder public: template - void operator()(VecField const& B, VecField const& E, VecField& Bnew, double dt) + void operator()(VecField const& B, VecField const& E, VecField& Bnew, floater_t<4> const dt) { if (!this->hasLayout()) throw std::runtime_error( @@ -49,7 +49,7 @@ class Faraday : public LayoutHolder private: - double dt_; + floater_t<4> dt_; template diff --git a/src/core/numerics/interpolator/interpolator.hpp b/src/core/numerics/interpolator/interpolator.hpp index a0a6a5a96..993f6e4a5 100644 --- a/src/core/numerics/interpolator/interpolator.hpp +++ b/src/core/numerics/interpolator/interpolator.hpp @@ -57,12 +57,15 @@ class Weighter template<> class Weighter<1> { + constexpr static floater_t<0> one = 1; + public: - inline void computeWeight(double normalizedPos, int startIndex, - std::array& weights) + template + inline void computeWeight(T normalizedPos, int startIndex, + std::array& weights) { - weights[1] = normalizedPos - static_cast(startIndex); - weights[0] = 1. - weights[1]; + weights[1] = normalizedPos - static_cast(startIndex); + weights[0] = one - weights[1]; } static constexpr int interp_order = 1; @@ -74,20 +77,24 @@ class Weighter<1> template<> class Weighter<2> { + constexpr static floater_t<0> p5 = .5; + constexpr static floater_t<0> p75 = .75; + public: - inline void computeWeight(double normalizedPos, int startIndex, - std::array& weights) + template + inline void computeWeight(T normalizedPos, int startIndex, + std::array& weights) { auto index = startIndex + 1; - auto delta = static_cast(index) - normalizedPos; - double coef1, coef2, coef3; - coef1 = 0.5 + delta; + auto delta = static_cast(index) - normalizedPos; + T coef1, coef2, coef3; + coef1 = p5 + delta; coef2 = delta; - coef3 = 0.5 - delta; + coef3 = p5 - delta; - weights[0] = 0.5 * coef1 * coef1; - weights[1] = 0.75 - coef2 * coef2; - weights[2] = 0.5 * coef3 * coef3; + weights[0] = p5 * coef1 * coef1; + weights[1] = p75 - coef2 * coef2; + weights[2] = p5 * coef3 * coef3; } static constexpr int interp_order = 2; @@ -100,27 +107,34 @@ class Weighter<2> template<> class Weighter<3> { + constexpr static floater_t<0> _1 = 1.0; + constexpr static floater_t<0> _2 = 2.0; + constexpr static floater_t<0> _3 = 3.0; + constexpr static floater_t<0> _4 = 4.0; + constexpr static floater_t<0> p5 = .5; + public: - inline void computeWeight(double normalizedPos, int startIndex, - std::array& weights) + template + inline void computeWeight(T normalizedPos, int startIndex, + std::array& weights) { - constexpr double _4_over_3 = 4. / 3.; - constexpr double _2_over_3 = 2. / 3.; + constexpr floater_t<0> _4_over_3 = _4 / _3; + constexpr floater_t<0> _2_over_3 = _2 / _3; - auto index = static_cast(startIndex) - normalizedPos; - double coef1 = 1. + 0.5 * index; - double coef2 = index + 1; - double coef3 = index + 2; - double coef4 = 1. - 0.5 * (index + 3); + auto index = static_cast>(startIndex) - normalizedPos; + floater_t<0> const coef1 = _1 + p5 * index; + floater_t<0> const coef2 = index + 1; + floater_t<0> const coef3 = index + 2; + floater_t<0> const coef4 = _1 - p5 * (index + 3); - double coef2_sq = coef2 * coef2; - double coef2_cub = coef2_sq * coef2; - double coef3_sq = coef3 * coef3; - double coef3_cub = coef3_sq * coef3; + floater_t<0> const coef2_sq = coef2 * coef2; + floater_t<0> const coef2_cub = coef2_sq * coef2; + floater_t<0> const coef3_sq = coef3 * coef3; + floater_t<0> const coef3_cub = coef3_sq * coef3; weights[0] = _4_over_3 * coef1 * coef1 * coef1; - weights[1] = _2_over_3 - coef2_sq - 0.5 * coef2_cub; - weights[2] = _2_over_3 - coef3_sq + 0.5 * coef3_cub; + weights[1] = _2_over_3 - coef2_sq - p5 * coef2_cub; + weights[2] = _2_over_3 - coef3_sq + p5 * coef3_cub; weights[3] = _4_over_3 * coef4 * coef4 * coef4; } @@ -168,8 +182,8 @@ class MeshToParticle<1> auto const& [xStartIndex, xWeights] = start_index_and_weights_for_qty<0, GridLayout, quantity>(indexWeights); - auto const& order_size = xWeights.size(); - auto fieldAtParticle = 0.; + auto const& order_size = xWeights.size(); + floater_t<4> fieldAtParticle = 0.; for (auto ik = 0u; ik < order_size; ++ik) { @@ -203,11 +217,11 @@ class MeshToParticle<2> auto const& order_size = xWeights.size(); - double fieldAtParticle = 0.; + floater_t<4> fieldAtParticle = 0.; for (auto ix = 0u; ix < order_size; ++ix) { - double Yinterp = 0.; + floater_t<4> Yinterp = 0.; for (auto iy = 0u; iy < order_size; ++iy) { Yinterp += field(xStartIndex + ix, yStartIndex + iy) * yWeights[iy]; @@ -246,13 +260,13 @@ class MeshToParticle<3> auto const& order_size = xWeights.size(); - double fieldAtParticle = 0.; + floater_t<4> fieldAtParticle = 0.; for (auto ix = 0u; ix < order_size; ++ix) { - double Yinterp = 0.; + floater_t<4> Yinterp = 0.; for (auto iy = 0u; iy < order_size; ++iy) { - double Zinterp = 0.; + floater_t<4> Zinterp = 0.; for (auto iz = 0u; iz < order_size; ++iz) { Zinterp += field(xStartIndex + ix, yStartIndex + iy, zStartIndex + iz) @@ -284,7 +298,8 @@ class ParticleToMesh<1> public: template inline void operator()(Field& field, Particle const& particle, Func&& func, - Indexes const& startIndex, Weights const& weights, double coef = 1.) + Indexes const& startIndex, Weights const& weights, + floater_t<0> const coef = 1.) { auto const& [xStartIndex] = startIndex; auto const& [xWeights] = weights; @@ -309,7 +324,8 @@ class ParticleToMesh<2> public: template inline void operator()(Field& field, Particle const& particle, Func&& func, - Indexes const& startIndex, Weights const& weights, double coef = 1.) + Indexes const& startIndex, Weights const& weights, + floater_t<0> const coef = 1.) { auto const& [xStartIndex, yStartIndex] = startIndex; auto const& [xWeights, yWeights] = weights; @@ -321,8 +337,8 @@ class ParticleToMesh<2> { for (auto iy = 0u; iy < order_size; ++iy) { - auto x = xStartIndex + ix; - auto y = yStartIndex + iy; + auto const x = xStartIndex + ix; + auto const y = yStartIndex + iy; field(x, y) += deposit * xWeights[ix] * yWeights[iy]; } @@ -340,7 +356,8 @@ class ParticleToMesh<3> public: template inline void operator()(Field& field, Particle const& particle, Func&& func, - Indexes const& startIndex, Weights const& weights, double coef = 1.) + Indexes const& startIndex, Weights const& weights, + floater_t<0> const coef = 1.) { auto const& [xStartIndex, yStartIndex, zStartIndex] = startIndex; auto const& [xWeights, yWeights, zWeights] = weights; @@ -354,9 +371,9 @@ class ParticleToMesh<3> { for (auto iz = 0u; iz < order_size; ++iz) { - auto x = xStartIndex + ix; - auto y = yStartIndex + iy; - auto z = zStartIndex + iz; + auto const x = xStartIndex + ix; + auto const y = yStartIndex + iy; + auto const z = zStartIndex + iz; field(x, y, z) += deposit * xWeights[ix] * yWeights[iy] * zWeights[iz]; } @@ -384,7 +401,7 @@ class Interpolator : private Weighter auto indexAndWeights_(GridLayout const& layout, ICell const& iCell_, Delta const& delta) { // dual weights require -.5 to take the correct position weight - auto constexpr dual_offset = .5; + [[maybe_unused]] floater_t<0> constexpr dual_offset = .5; auto const& [startIndex_, weights_] = [&]() { if constexpr (centering == QtyCentering::dual) @@ -399,7 +416,7 @@ class Interpolator : private Weighter startIndex_[iDim] = iCell[iDim] - computeStartLeftShift(delta[iDim]); - double normalizedPos = iCell[iDim] + delta[iDim]; + floater_t<0> normalizedPos = iCell[iDim] + delta[iDim]; if constexpr (centering == QtyCentering::dual) normalizedPos -= dual_offset; @@ -409,6 +426,7 @@ class Interpolator : private Weighter } public: + floater_t<0> static constexpr one = 1; auto static constexpr interp_order = interpOrder; auto static constexpr dimension = dim; @@ -423,7 +441,7 @@ class Interpolator : private Weighter template inline auto operator()(Particle_t& currPart, Electromag const& Em, GridLayout const& layout) { - using E_B_tuple = std::tuple, std::array>; + using E_B_tuple = std::tuple, 3>, std::array, 3>>; using Scalar = HybridQuantity::Scalar; // for each particle, first calculate the startIndex and weights for dual and @@ -470,7 +488,7 @@ class Interpolator : private Weighter */ template inline void operator()(ParticleRange& particleRange, Field& density, VecField& flux, - GridLayout const& layout, double coef = 1.) + GridLayout const& layout, floater_t<0> const coef = 1.) { auto begin = particleRange.begin(); auto end = particleRange.end(); @@ -487,7 +505,7 @@ class Interpolator : private Weighter currPart->delta); particleToMesh_( - density, *currPart, [](auto const& part) { return 1.; }, startIndex_, weights_, + density, *currPart, [](auto const& part) { return one; }, startIndex_, weights_, coef); particleToMesh_( xFlux, *currPart, [](auto const& part) { return part.v[0]; }, startIndex_, weights_, @@ -503,7 +521,7 @@ class Interpolator : private Weighter } template inline void operator()(ParticleRange&& range, Field& density, VecField& flux, - GridLayout const& layout, double coef = 1.) + GridLayout const& layout, floater_t<0> const coef = 1.) { (*this)(range, density, flux, layout, coef); } @@ -552,7 +570,7 @@ class Interpolator : private Weighter static_assert(dimension <= 3 && dimension > 0 && interpOrder >= 1 && interpOrder <= 3, "error"); using Starts = std::array; - using Weights = std::array, dimension>; + using Weights = std::array, nbrPointsSupport(interpOrder)>, dimension>; Weighter weightComputer_; MeshToParticle meshToParticle_; diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 27d265e6d..8ac088c2c 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -22,8 +22,8 @@ class Ohm : public LayoutHolder public: explicit Ohm(PHARE::initializer::PHAREDict const& dict) - : eta_{dict["resistivity"].template to()} - , nu_{dict["hyper_resistivity"].template to()} + : eta_{dict["resistivity"].template to>()} + , nu_{dict["hyper_resistivity"].template to>()} , hyper_mode{cppdict::get_value(dict, "hyper_mode", std::string{"constant"}) == "constant" ? HyperMode::constant : HyperMode::spatial} @@ -57,8 +57,8 @@ class Ohm : public LayoutHolder private: - double const eta_; - double const nu_; + floater_t<4> const eta_; + floater_t<4> const nu_; HyperMode const hyper_mode; template @@ -87,7 +87,7 @@ class Ohm : public LayoutHolder template - auto ideal1D_(VecField const& Ve, VecField const& B, MeshIndex<1> index) const + floater_t<4> ideal1D_(VecField const& Ve, VecField const& B, MeshIndex<1> index) const { if constexpr (component == Component::X) { @@ -250,7 +250,7 @@ class Ohm : public LayoutHolder template - auto ideal_(VecField const& Ve, VecField const& B, MeshIndex index) const + floater_t<4> ideal_(VecField const& Ve, VecField const& B, MeshIndex index) const { if constexpr (dimension == 1) return ideal1D_(Ve, B, index); @@ -360,17 +360,20 @@ class Ohm : public LayoutHolder auto spatial_hyperresistive_(VecField const& J, VecField const& B, Field const& n, MeshIndex index) const { // TODO : https://github.com/PHAREHUB/PHARE/issues/3 - auto const lvlCoeff = 1. / std::pow(4, layout_->levelNumber()); - auto constexpr min_density = 0.1; + auto const lvlCoeff + = 1.f / static_cast>(std::pow(4.f, layout_->levelNumber())); + auto constexpr min_density = 0.1f; auto computeHR = [&](auto BxProj, auto ByProj, auto BzProj, auto nProj) { - auto const BxOnE = GridLayout::project(B(Component::X), index, BxProj); - auto const ByOnE = GridLayout::project(B(Component::Y), index, ByProj); - auto const BzOnE = GridLayout::project(B(Component::Z), index, BzProj); - auto const nOnE = GridLayout::project(n, index, nProj); - auto b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE); - return -nu_ * (b / (nOnE + min_density) + 1) * lvlCoeff - * layout_->laplacian(J(component), index); + auto const BxOnE = GridLayout::project(B(Component::X), index, BxProj); + auto const ByOnE = GridLayout::project(B(Component::Y), index, ByProj); + auto const BzOnE = GridLayout::project(B(Component::Z), index, BzProj); + auto const nOnE = GridLayout::project(n, index, nProj); + floater_t<4> const b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE); + floater_t<4> const ret = -nu_ * (b / (nOnE + min_density) + 1) * lvlCoeff + * layout_->laplacian(J(component), index); + static_assert(std::is_same_v, floater_t<4>>); + return ret; }; if constexpr (component == Component::X) diff --git a/src/core/numerics/pusher/boris.hpp b/src/core/numerics/pusher/boris.hpp index 0d8ff09e8..46fd30d4f 100644 --- a/src/core/numerics/pusher/boris.hpp +++ b/src/core/numerics/pusher/boris.hpp @@ -1,17 +1,17 @@ #ifndef PHARE_CORE_PUSHER_BORIS_HPP #define PHARE_CORE_PUSHER_BORIS_HPP +#include "core/errors.hpp" +#include "core/logger.hpp" +#include "core/numerics/pusher/pusher.hpp" + + #include #include #include -#include #include -#include -#include "core/numerics/pusher/pusher.hpp" -#include "core/utilities/range/range.hpp" -#include "core/errors.hpp" -#include "core/logger.hpp" -#include "core/data/particles/particle.hpp" +#include + namespace PHARE::core { @@ -35,7 +35,7 @@ class BorisPusher /** see Pusher::move() documentation*/ #if 0 ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, Interpolator& interpolator, + Electromag const& emFields, floater_t<4> mass, Interpolator& interpolator, ParticleSelector const& particleIsNotLeaving, BoundaryCondition& bc, GridLayout const& layout) override { @@ -77,7 +77,7 @@ class BorisPusher ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, Interpolator& interpolator, + Electromag const& emFields, floater_t<4> mass, Interpolator& interpolator, GridLayout const& layout, ParticleSelector firstSelector, ParticleSelector secondSelector) override { @@ -92,7 +92,7 @@ class BorisPusher rangeOut = firstSelector(rangeOut); - double const dto2m = 0.5 * dt_ / mass; + floater_t<4> const dto2m = 0.5f * dt_ / mass; for (auto idx = rangeOut.ibegin(); idx < rangeOut.iend(); ++idx) { auto& currPart = rangeOut.array()[idx]; @@ -112,10 +112,11 @@ class BorisPusher /** see Pusher::move() documentation*/ - void setMeshAndTimeStep(std::array ms, double const ts) override + void setMeshAndTimeStep(std::array, dim> const ms, floater_t<4> const ts) override { + floater_t<4> constexpr static p5 = .5; std::transform(std::begin(ms), std::end(ms), std::begin(halfDtOverDl_), - [ts](double& x) { return 0.5 * ts / x; }); + [ts](auto& x) { return p5 * ts / x; }); dt_ = ts; } @@ -127,13 +128,14 @@ class BorisPusher template auto advancePosition_(Particle const& partIn, Particle& partOut) { + using Float_t = floater_t<0>; + std::array newCell; for (std::size_t iDim = 0; iDim < dim; ++iDim) { - double delta - = partIn.delta[iDim] + static_cast(halfDtOverDl_[iDim] * partIn.v[iDim]); + Float_t delta = partIn.delta[iDim] + halfDtOverDl_[iDim] * Float_t(partIn.v[iDim]); - double iCell = std::floor(delta); + Float_t iCell = std::floor(delta); if (std::abs(delta) > 2) { PHARE_LOG_ERROR("Error, particle moves more than 1 cell, delta >2"); @@ -190,55 +192,55 @@ class BorisPusher /** Accelerate the particles in rangeIn and put the new velocity in rangeOut */ template - void accelerate_(Particle_t& part, ParticleEB const& particleEB, double const& dto2m) + void accelerate_(Particle_t& part, ParticleEB const& particleEB, floater_t<4> const& dto2m) { auto& [pE, pB] = particleEB; auto& [pEx, pEy, pEz] = pE; auto& [pBx, pBy, pBz] = pB; - double const coef1 = part.charge * dto2m; + floater_t<4> const coef1 = part.charge * dto2m; // We now apply the 3 steps of the BORIS PUSHER // 1st half push of the electric field - double velx1 = part.v[0] + coef1 * pEx; - double vely1 = part.v[1] + coef1 * pEy; - double velz1 = part.v[2] + coef1 * pEz; + floater_t<4> velx1 = part.v[0] + coef1 * pEx; + floater_t<4> vely1 = part.v[1] + coef1 * pEy; + floater_t<4> velz1 = part.v[2] + coef1 * pEz; // preparing variables for magnetic rotation - double const rx = coef1 * pBx; - double const ry = coef1 * pBy; - double const rz = coef1 * pBz; + floater_t<4> const rx = coef1 * pBx; + floater_t<4> const ry = coef1 * pBy; + floater_t<4> const rz = coef1 * pBz; - double const rx2 = rx * rx; - double const ry2 = ry * ry; - double const rz2 = rz * rz; - double const rxry = rx * ry; - double const rxrz = rx * rz; - double const ryrz = ry * rz; + floater_t<4> const rx2 = rx * rx; + floater_t<4> const ry2 = ry * ry; + floater_t<4> const rz2 = rz * rz; + floater_t<4> const rxry = rx * ry; + floater_t<4> const rxrz = rx * rz; + floater_t<4> const ryrz = ry * rz; - double const invDet = 1. / (1. + rx2 + ry2 + rz2); + floater_t<4> const invDet = 1.f / (1.f + rx2 + ry2 + rz2); // preparing rotation matrix due to the magnetic field // m = invDet*(I + r*r - r x I) - I where x denotes the cross product - double const mxx = 1. + rx2 - ry2 - rz2; - double const mxy = 2. * (rxry + rz); - double const mxz = 2. * (rxrz - ry); + floater_t<4> const mxx = 1.f + rx2 - ry2 - rz2; + floater_t<4> const mxy = 2.f * (rxry + rz); + floater_t<4> const mxz = 2.f * (rxrz - ry); - double const myx = 2. * (rxry - rz); - double const myy = 1. + ry2 - rx2 - rz2; - double const myz = 2. * (ryrz + rx); + floater_t<4> const myx = 2.f * (rxry - rz); + floater_t<4> const myy = 1.f + ry2 - rx2 - rz2; + floater_t<4> const myz = 2.f * (ryrz + rx); - double const mzx = 2. * (rxrz + ry); - double const mzy = 2. * (ryrz - rx); - double const mzz = 1. + rz2 - rx2 - ry2; + floater_t<4> const mzx = 2.f * (rxrz + ry); + floater_t<4> const mzy = 2.f * (ryrz - rx); + floater_t<4> const mzz = 1.f + rz2 - rx2 - ry2; // magnetic rotation - double const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet; - double const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet; - double const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet; + floater_t<4> const velx2 = (mxx * velx1 + mxy * vely1 + mxz * velz1) * invDet; + floater_t<4> const vely2 = (myx * velx1 + myy * vely1 + myz * velz1) * invDet; + floater_t<4> const velz2 = (mzx * velx1 + mzy * vely1 + mzz * velz1) * invDet; // 2nd half push of the electric field @@ -255,8 +257,8 @@ class BorisPusher - std::array halfDtOverDl_; - double dt_; + std::array, dim> halfDtOverDl_; + floater_t<4> dt_; }; } // namespace PHARE::core diff --git a/src/core/numerics/pusher/pusher.hpp b/src/core/numerics/pusher/pusher.hpp index 9c1b48b9f..1ecf3fea8 100644 --- a/src/core/numerics/pusher/pusher.hpp +++ b/src/core/numerics/pusher/pusher.hpp @@ -1,14 +1,10 @@ #ifndef PHARE_CORE_NUMERICS_PUSHER_PUSHER_HPP #define PHARE_CORE_NUMERICS_PUSHER_PUSHER_HPP + #include -#include -#include #include -#include "core/utilities/range/range.hpp" -#include "core/data/particles/particle.hpp" - namespace PHARE { namespace core @@ -26,13 +22,15 @@ namespace core // TODO : to really be independant on boris which has 2 push steps // we should have an arbitrary number of selectors, 1 per push step virtual ParticleRange move(ParticleRange const& rangeIn, ParticleRange& rangeOut, - Electromag const& emFields, double mass, + Electromag const& emFields, floater_t<4> mass, Interpolator& interpolator, GridLayout const& layout, ParticleSelector firstSelector, ParticleSelector secondSelector) = 0; - virtual void setMeshAndTimeStep(std::array ms, double ts) = 0; + virtual void setMeshAndTimeStep(std::array, dim> const ms, + floater_t<4> const ts) + = 0; virtual ~Pusher() {} }; diff --git a/src/core/utilities/logger/logger_defaults.hpp b/src/core/utilities/logger/logger_defaults.hpp index 6bb2b4c90..5ae689122 100644 --- a/src/core/utilities/logger/logger_defaults.hpp +++ b/src/core/utilities/logger/logger_defaults.hpp @@ -1,3 +1,5 @@ +// IWYU pragma: private, include "core/logger.hpp" + #ifndef PHARE_CORE_UTILITIES_LOGGER_LOGGER_DEFAULTS_HPP #define PHARE_CORE_UTILITIES_LOGGER_LOGGER_DEFAULTS_HPP diff --git a/src/diagnostic/diagnostic_model_view.hpp b/src/diagnostic/diagnostic_model_view.hpp index 1dafc2826..4e1524404 100644 --- a/src/diagnostic/diagnostic_model_view.hpp +++ b/src/diagnostic/diagnostic_model_view.hpp @@ -37,8 +37,8 @@ class ModelView : public IModelView using GridLayout = typename Model::gridlayout_type; using PatchProperties = cppdict::Dict, std::vector, - std::vector, std::vector, std::string, - std::vector>; + std::vector, std::vector, std::vector, + std::string, std::vector>; ModelView(Hierarchy& hierarchy, Model& model) : model_{model} diff --git a/src/initializer/data_provider.hpp b/src/initializer/data_provider.hpp index c48f31a6e..db73595db 100644 --- a/src/initializer/data_provider.hpp +++ b/src/initializer/data_provider.hpp @@ -21,36 +21,36 @@ namespace initializer }; template<> - struct InitFunctionHelper + struct InitFunctionHelper, 1> { - using return_type = std::shared_ptr>; - using param_type = std::vector const&; + using return_type = std::shared_ptr>>; + using param_type = std::vector> const&; using type = std::function; }; template<> - struct InitFunctionHelper + struct InitFunctionHelper, 2> { - using return_type = std::shared_ptr>; - using param_type = std::vector const&; + using return_type = std::shared_ptr>>; + using param_type = std::vector> const&; using type = std::function; }; template<> - struct InitFunctionHelper + struct InitFunctionHelper, 3> { - using return_type = std::shared_ptr>; - using param_type = std::vector const&; + using return_type = std::shared_ptr>>; + using param_type = std::vector> const&; using type = std::function; }; template - using InitFunction = typename InitFunctionHelper::type; + using InitFunction = typename InitFunctionHelper, dim>::type; - using PHAREDict = cppdict::Dict, double, std::vector, - std::size_t, std::optional, std::string, - InitFunction<1>, InitFunction<2>, InitFunction<3>>; + using PHAREDict = cppdict::Dict, float, double, std::vector, + std::vector, std::size_t, std::optional, + std::string, InitFunction<1>, InitFunction<2>, InitFunction<3>>; class PHAREDictHandler diff --git a/src/initializer/dictator.cpp b/src/initializer/dictator.cpp index 5d41fa7e8..73fc49314 100644 --- a/src/initializer/dictator.cpp +++ b/src/initializer/dictator.cpp @@ -49,12 +49,14 @@ PYBIND11_MODULE(dictator, m) m.def("add_int", add, "add"); m.def("add_vector_int", add>, "add"); m.def("add_double", add, "add"); + m.def("add_double", add, "add"); m.def("add_string", add, "add"); m.def("addInitFunction1D", add>, "add"); m.def("addInitFunction2D", add>, "add"); m.def("addInitFunction3D", add>, "add"); + m.def("add_array_as_vector", add_array_as_vector, "add_array_as_vector"); m.def("add_array_as_vector", add_array_as_vector, "add_array_as_vector"); m.def("stop", []() { PHARE::initializer::PHAREDictHandler::INSTANCE().stop(); }); diff --git a/src/phare_core.hpp b/src/phare_core.hpp index 0d6ef0d83..536e24483 100644 --- a/src/phare_core.hpp +++ b/src/phare_core.hpp @@ -1,30 +1,23 @@ #ifndef PHARE_CORE_INCLUDE_HPP #define PHARE_CORE_INCLUDE_HPP +#include "core/data/ions/ions.hpp" #include "core/data/grid/grid.hpp" -#include "core/data/electromag/electromag.hpp" -#include "core/data/electrons/electrons.hpp" #include "core/data/grid/gridlayout.hpp" +#include "core/data/vecfield/vecfield.hpp" +#include "core/data/electrons/electrons.hpp" +#include "core/data/electromag/electromag.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" #include "core/data/grid/gridlayoutimplyee.hpp" +#include "core/data/particles/particle_array.hpp" #include "core/data/ions/ion_population/ion_population.hpp" -#include "core/data/ions/ions.hpp" #include "core/data/ions/particle_initializers/maxwellian_particle_initializer.hpp" -#include "core/data/ndarray/ndarray_vector.hpp" -#include "core/data/particles/particle_array.hpp" -#include "core/data/vecfield/vecfield.hpp" -#include "core/models/physical_state.hpp" -#include "core/models/physical_state.hpp" -#include "core/utilities/meta/meta_utilities.hpp" -#include "core/utilities/algorithm.hpp" -#include "core/logger.hpp" -#include -#include -#include +#include "cppdict/include/dict.hpp" + #include #include -#include "cppdict/include/dict.hpp" namespace PHARE::core { @@ -34,11 +27,11 @@ struct PHARE_Types static auto constexpr dimension = dimension_; static auto constexpr interp_order = interp_order_; - using Array_t = PHARE::core::NdArrayVector; - using ArrayView_t = PHARE::core::NdArrayView; - using Grid_t = PHARE::core::Grid; - using Field_t = PHARE::core::Field; - using VecField_t = PHARE::core::VecField; + using Array_t = PHARE::core::NdArrayVector>; + using ArrayView_t = PHARE::core::NdArrayView>; + using Grid_t = PHARE::core::Grid; + using Field_t = core::Field>; + using VecField_t = PHARE::core::VecField; using SymTensorField_t = PHARE::core::SymTensorField; using Electromag_t = PHARE::core::Electromag; using YeeLayout_t = PHARE::core::GridLayoutImplYee; diff --git a/src/python3/cpp_simulator.cpp b/src/python3/cpp_simulator.cpp index a5c0a5d1c..2d9615e11 100644 --- a/src/python3/cpp_simulator.cpp +++ b/src/python3/cpp_simulator.cpp @@ -19,8 +19,8 @@ PYBIND11_MODULE(PHARE_CPP_MOD_NAME, m) core::apply(core::possibleSimulators(), [&](auto const& simType) { declare_all(m, simType); }); - declarePatchData, 1>(m, "PatchDataVectorDouble_1D"); - declarePatchData, 2>(m, "PatchDataVectorDouble_2D"); - declarePatchData, 3>(m, "PatchDataVectorDouble_3D"); + // declarePatchData>, 1>(m, "PatchDataVectorDouble_1D"); + // declarePatchData>, 2>(m, "PatchDataVectorDouble_2D"); + // declarePatchData>, 3>(m, "PatchDataVectorDouble_3D"); } } // namespace PHARE::pydata diff --git a/src/python3/cpp_simulator.hpp b/src/python3/cpp_simulator.hpp index 7b78094bb..4837640e3 100644 --- a/src/python3/cpp_simulator.hpp +++ b/src/python3/cpp_simulator.hpp @@ -21,9 +21,9 @@ #include "pybind11/functional.h" #include "python3/particles.hpp" -#include "python3/patch_data.hpp" -#include "python3/patch_level.hpp" -#include "python3/data_wrangler.hpp" +// #include "python3/patch_data.hpp" +// #include "python3/patch_level.hpp" +// #include "python3/data_wrangler.hpp" @@ -31,18 +31,18 @@ namespace py = pybind11; namespace PHARE::pydata { -template -void declarePatchData(py::module& m, std::string key) -{ - using PatchDataType = PatchData; - py::class_(m, key.c_str()) - .def_readonly("patchID", &PatchDataType::patchID) - .def_readonly("origin", &PatchDataType::origin) - .def_readonly("lower", &PatchDataType::lower) - .def_readonly("upper", &PatchDataType::upper) - .def_readonly("nGhosts", &PatchDataType::nGhosts) - .def_readonly("data", &PatchDataType::data); -} +// template +// void declarePatchData(py::module& m, std::string key) +// { +// using PatchDataType = PatchData; +// py::class_(m, key.c_str()) +// .def_readonly("patchID", &PatchDataType::patchID) +// .def_readonly("origin", &PatchDataType::origin) +// .def_readonly("lower", &PatchDataType::lower) +// .def_readonly("upper", &PatchDataType::upper) +// .def_readonly("nGhosts", &PatchDataType::nGhosts) +// .def_readonly("data", &PatchDataType::data); +// } template void declareDim(py::module& m) @@ -58,8 +58,8 @@ void declareDim(py::module& m) .def_readwrite("v", &CP::v) .def("size", &CP::size); - name = "PatchData" + name; - declarePatchData(m, name.c_str()); + // name = "PatchData" + name; + // declarePatchData(m, name.c_str()); } template @@ -87,40 +87,41 @@ void declare_etc(py::module& m) std::string type_string = "_" + std::to_string(dim) + "_" + std::to_string(interp) + "_" + std::to_string(nbRefinedPart); - using Sim = Simulator; - using DW = DataWrangler; - std::string name = "DataWrangler" + type_string; - py::class_>(m, name.c_str()) - .def(py::init const&, std::shared_ptr const&>()) - .def(py::init const&, std::shared_ptr const&>()) - .def("sync_merge", &DW::sync_merge) - .def("getPatchLevel", &DW::getPatchLevel) - .def("getNumberOfLevels", &DW::getNumberOfLevels); - - using PL = PatchLevel; - name = "PatchLevel_" + type_string; - - py::class_>(m, name.c_str()) - .def("getEM", &PL::getEM) - .def("getE", &PL::getE) - .def("getB", &PL::getB) - .def("getBx", &PL::getBx) - .def("getBy", &PL::getBy) - .def("getBz", &PL::getBz) - .def("getEx", &PL::getEx) - .def("getEy", &PL::getEy) - .def("getEz", &PL::getEz) - .def("getVix", &PL::getVix) - .def("getViy", &PL::getViy) - .def("getViz", &PL::getViz) - .def("getDensity", &PL::getDensity) - .def("getBulkVelocity", &PL::getBulkVelocity) - .def("getPopDensities", &PL::getPopDensities) - .def("getPopFluxes", &PL::getPopFlux) - .def("getFx", &PL::getFx) - .def("getFy", &PL::getFy) - .def("getFz", &PL::getFz) - .def("getParticles", &PL::getParticles, py::arg("userPopName") = "all"); + // using Sim = Simulator; + + std::string name; + + // using DW = DataWrangler; + // py::class_>(m, name.c_str()) + // .def(py::init const&, std::shared_ptr const&>()) + // .def(py::init const&, std::shared_ptr + // const&>()) .def("sync_merge", &DW::sync_merge) .def("getPatchLevel", &DW::getPatchLevel) + // .def("getNumberOfLevels", &DW::getNumberOfLevels); + + // using PL = PatchLevel; + // name = "PatchLevel_" + type_string; + + // py::class_>(m, name.c_str()) + // .def("getEM", &PL::getEM) + // .def("getE", &PL::getE) + // .def("getB", &PL::getB) + // .def("getBx", &PL::getBx) + // .def("getBy", &PL::getBy) + // .def("getBz", &PL::getBz) + // .def("getEx", &PL::getEx) + // .def("getEy", &PL::getEy) + // .def("getEz", &PL::getEz) + // .def("getVix", &PL::getVix) + // .def("getViy", &PL::getViy) + // .def("getViz", &PL::getViz) + // .def("getDensity", &PL::getDensity) + // .def("getBulkVelocity", &PL::getBulkVelocity) + // .def("getPopDensities", &PL::getPopDensities) + // .def("getPopFluxes", &PL::getPopFlux) + // .def("getFx", &PL::getFx) + // .def("getFy", &PL::getFy) + // .def("getFz", &PL::getFz) + // .def("getParticles", &PL::getParticles, py::arg("userPopName") = "all"); using _Splitter = PHARE::amr::Splitter<_dim, _interp, core::RefinedParticlesConst>; diff --git a/src/python3/particles.hpp b/src/python3/particles.hpp index 84d0e2411..f582f14aa 100644 --- a/src/python3/particles.hpp +++ b/src/python3/particles.hpp @@ -16,21 +16,21 @@ namespace PHARE::pydata template core::ContiguousParticlesView contiguousViewFrom(PyArrayTuple const& py_particles) { - return {makeSpan(std::get<0>(py_particles)), // iCell - makeSpan(std::get<1>(py_particles)), // delta - makeSpan(std::get<2>(py_particles)), // weight - makeSpan(std::get<3>(py_particles)), // charge - makeSpan(std::get<4>(py_particles))}; // v + return {makeSpan(std::get<0>(py_particles)), // iCell + makeSpan>(std::get<1>(py_particles)), // delta + makeSpan>(std::get<2>(py_particles)), // weight + makeSpan>(std::get<3>(py_particles)), // charge + makeSpan>(std::get<4>(py_particles))}; // v } template pyarray_particles_t makePyArrayTuple(std::size_t const size) { - return std::make_tuple(py_array_t(size * dim), // iCell - py_array_t(size * dim), // delta - py_array_t(size), // weight - py_array_t(size), // charge - py_array_t(size * 3)); // v + return std::make_tuple(py_array_t(size * dim), // iCell + py_array_t>(size * dim), // delta + py_array_t>(size), // weight + py_array_t>(size), // charge + py_array_t>(size * 3)); // v } diff --git a/src/python3/pybind_def.hpp b/src/python3/pybind_def.hpp index 256b5d6c9..8265bbbc4 100644 --- a/src/python3/pybind_def.hpp +++ b/src/python3/pybind_def.hpp @@ -1,12 +1,15 @@ #ifndef PHARE_PYTHON_PYBIND_DEF_HPP #define PHARE_PYTHON_PYBIND_DEF_HPP +#include "core/def.hpp" +#include "core/utilities/span.hpp" + #include #include #include #include -#include "core/utilities/span.hpp" + #include "pybind11/stl.h" #include "pybind11/numpy.h" @@ -19,12 +22,19 @@ template using py_array_t = pybind11::array_t; -using pyarray_particles_t = std::tuple, py_array_t, py_array_t, - py_array_t, py_array_t>; +using pyarray_particles_t = std::tuple, // iCell + py_array_t>, // delta + py_array_t>, // weight + py_array_t>, // charge + py_array_t> // v + >; -using pyarray_particles_crt - = std::tuple const&, py_array_t const&, py_array_t const&, - py_array_t const&, py_array_t const&>; +using pyarray_particles_crt = std::tuple const&, // iCell + py_array_t> const&, // delta + py_array_t> const&, // weight + py_array_t> const&, // charge + py_array_t> const& // v + >; template std::size_t ndSize(PyArrayInfo const& ar_info) diff --git a/tests/amr/data/field/coarsening/test_linear_coarsen.hpp b/tests/amr/data/field/coarsening/test_linear_coarsen.hpp index c2a45c037..414d96d1c 100644 --- a/tests/amr/data/field/coarsening/test_linear_coarsen.hpp +++ b/tests/amr/data/field/coarsening/test_linear_coarsen.hpp @@ -1,13 +1,16 @@ #ifndef PHARE_TEST_LINEAR_COARSEN_HPP #define PHARE_TEST_LINEAR_COARSEN_HPP +#include "phare_core.hpp" + #include "amr/data/field/coarsening/default_field_coarsener.hpp" #include "amr/data/field/coarsening/magnetic_field_coarsener.hpp" #include "amr/data/field/coarsening/field_coarsen_operator.hpp" #include "amr/data/field/coarsening/field_coarsen_index_weight.hpp" -#include "core/data/grid/grid.hpp" -#include "core/data/grid/gridlayout.hpp" -#include "core/data/grid/gridlayout_impl.hpp" + +// #include "core/data/grid/grid.hpp" +// #include "core/data/grid/gridlayout.hpp" +// #include "core/data/grid/gridlayout_impl.hpp" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -52,11 +55,12 @@ struct Files } }; -template +template struct EMData { - using Grid_t = Grid, HybridQuantity::Scalar>; - using GridPtr = std::shared_ptr; + using PHARE_Types = PHARE::core::PHARE_Types; + using Grid_t = PHARE_Types::Grid_t; + using GridPtr = std::shared_ptr; std::string em_key; @@ -83,20 +87,22 @@ struct EMData * We get the fine and coarse data and the expected after coarse field * for two centering : primal and dual */ -template +template struct FieldCoarsenTestData { - static constexpr auto dimension = dimension_; + static constexpr auto dimension = dim; static constexpr auto interp = 1; static constexpr double absError = 1.e-8; - using GridYee_t = GridLayout>; - using Grid_t = Grid, HybridQuantity::Scalar>; + using PHARE_Types = PHARE::core::PHARE_Types; + using GridYee_t = PHARE_Types::GridLayout_t; + using Grid_t = PHARE_Types::Grid_t; + using Real_t = floater_t<4>; EMData em; - std::array meshSizeCoarse{ConstArray(0.2)}; - std::array meshSizeFine{ConstArray(0.1)}; + std::array meshSizeCoarse{ConstArray(0.2)}; + std::array meshSizeFine{ConstArray(0.1)}; std::array coarseIndexesX{{0, 39}}; std::array fineIndexesX{18, 37}; @@ -329,7 +335,10 @@ TYPED_TEST(FieldCoarsenOperatorTest, doTheExpectedCoarseningForEB) if constexpr (dim == 1) { for (auto ix = iStartX; ix <= iEndX; ++ix) - EXPECT_THAT(coarseValue(ix), DoubleNear(expectedCoarseValue(ix), absError)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(coarseValue(ix), DoubleNear(expectedCoarseValue(ix), absError)); + } } else { @@ -340,8 +349,11 @@ TYPED_TEST(FieldCoarsenOperatorTest, doTheExpectedCoarseningForEB) { for (auto ix = iStartX; ix <= iEndX; ++ix) for (auto iy = iStartY; iy <= iEndY; ++iy) - EXPECT_THAT(coarseValue(ix, iy), - DoubleNear(expectedCoarseValue(ix, iy), absError)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(coarseValue(ix, iy), + DoubleNear(expectedCoarseValue(ix, iy), absError)); + } } else if constexpr (dim == 3) // TODO 3D { diff --git a/tests/amr/data/field/copy_pack/field_data_test_param.hpp b/tests/amr/data/field/copy_pack/field_data_test_param.hpp index ab5b64d92..7df6b0e3c 100644 --- a/tests/amr/data/field/copy_pack/field_data_test_param.hpp +++ b/tests/amr/data/field/copy_pack/field_data_test_param.hpp @@ -278,7 +278,7 @@ struct AFieldData1DCenteredOnEy : public ::testing::Test // Using used later in test -using Grid1D = Grid, HybridQuantity::Scalar>; +using Grid1D = Grid>, HybridQuantity::Scalar>; using FieldDataTest1DOrder1 = FieldDataTestParam>, Grid1D>; using FieldDataTest1DOrder2 = FieldDataTestParam>, Grid1D>; diff --git a/tests/amr/data/field/geometry/test_field_geometry.cpp b/tests/amr/data/field/geometry/test_field_geometry.cpp index 0bb2ce95a..ca3c4dc45 100644 --- a/tests/amr/data/field/geometry/test_field_geometry.cpp +++ b/tests/amr/data/field/geometry/test_field_geometry.cpp @@ -16,6 +16,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include using testing::Eq; @@ -32,7 +33,7 @@ using namespace PHARE::amr; -using Grid1D = Grid, HybridQuantity::Scalar>; +using Grid1D = Grid>, HybridQuantity::Scalar>; template struct FieldGeometryParam diff --git a/tests/amr/data/field/refine/test_field_refine.cpp b/tests/amr/data/field/refine/test_field_refine.cpp index 836616d56..a2d25c743 100644 --- a/tests/amr/data/field/refine/test_field_refine.cpp +++ b/tests/amr/data/field/refine/test_field_refine.cpp @@ -78,7 +78,7 @@ TYPED_TEST(aFieldRefineOperator, canBeCreated) static constexpr auto interp = typename TypeParam::second_type{}(); using GridYee = GridLayout>; - using GridT = Grid, HybridQuantity::Scalar>; + using GridT = Grid>, HybridQuantity::Scalar>; FieldRefineOperator> linearRefine{}; } diff --git a/tests/amr/data/field/refine/test_field_refinement_on_hierarchy.cpp b/tests/amr/data/field/refine/test_field_refinement_on_hierarchy.cpp index d40e1c777..3dffb2abc 100644 --- a/tests/amr/data/field/refine/test_field_refinement_on_hierarchy.cpp +++ b/tests/amr/data/field/refine/test_field_refinement_on_hierarchy.cpp @@ -20,7 +20,7 @@ struct ALinearFieldRefineTest : public ::testing::Test static constexpr auto refine = 2; using GridYee = GridLayout>; - using GridND = Grid, HybridQuantity::Scalar>; + using GridND = Grid>, HybridQuantity::Scalar>; public: void SetUp() override diff --git a/tests/amr/data/field/time_interpolate/test_field_data_time_interpolate.cpp b/tests/amr/data/field/time_interpolate/test_field_data_time_interpolate.cpp index 80861eb61..657f87dd9 100644 --- a/tests/amr/data/field/time_interpolate/test_field_data_time_interpolate.cpp +++ b/tests/amr/data/field/time_interpolate/test_field_data_time_interpolate.cpp @@ -48,7 +48,7 @@ struct aFieldLinearTimeInterpolate : public ::testing::Test static constexpr auto interp = typename TypeInfo::second_type{}(); using GridYee = GridLayout>; - using GridND = Grid, HybridQuantity::Scalar>; + using GridND = Grid>, HybridQuantity::Scalar>; using FieldDataT = FieldData; FieldLinearTimeInterpolate timeOp{}; @@ -66,24 +66,26 @@ struct aFieldLinearTimeInterpolate : public ::testing::Test std::string const fieldName{"Bx"}; SAMRAI::hier::IntVector ghost{dimension, 5}; - static auto constexpr dl = ConstArray(0.01); + static auto constexpr dl = ConstArray, dim>(0.01); static constexpr auto nbrCells = ConstArray(upper - lower + 1); - Point origin{{0.}}; + Point, dim> const origin{{0.}}; + GridYee const layout; std::shared_ptr srcOld; std::shared_ptr srcNew; std::shared_ptr destNew; aFieldLinearTimeInterpolate() - : srcOld{std::make_shared(domain, ghost, fieldName, dl, nbrCells, origin, qty)} - , srcNew{std::make_shared(domain, ghost, fieldName, dl, nbrCells, origin, qty)} - , destNew{std::make_shared(domain, ghost, fieldName, dl, nbrCells, origin, qty)} + : layout{dl, nbrCells, origin} + , srcOld{std::make_shared(domain, ghost, fieldName, layout, qty)} + , srcNew{std::make_shared(domain, ghost, fieldName, layout, qty)} + , destNew{std::make_shared(domain, ghost, fieldName, layout, qty)} { - double oldTime = 0.; - double newTime = 0.5; + floater_t<4> oldTime = 0.; + floater_t<4> newTime = 0.5; srcOld->setTime(oldTime); srcNew->setTime(newTime); @@ -178,7 +180,7 @@ int aFieldLinearTimeInterpolate::countLocal = 0; TYPED_TEST(aFieldLinearTimeInterpolate, giveOldSrcForAlphaZero) { - double interpolateTime = 0.; + floater_t<4> interpolateTime = 0.; this->destNew->setTime(interpolateTime); auto& layout = this->srcOld->gridLayout; @@ -259,7 +261,7 @@ TYPED_TEST(aFieldLinearTimeInterpolate, giveOldSrcForAlphaZero) TYPED_TEST(aFieldLinearTimeInterpolate, giveNewSrcForAlphaOne) { - double interpolateTime = 0.5; + floater_t<4> interpolateTime = 0.5; this->destNew->setTime(interpolateTime); auto& layout = this->srcOld->gridLayout; @@ -339,7 +341,7 @@ TYPED_TEST(aFieldLinearTimeInterpolate, giveNewSrcForAlphaOne) TYPED_TEST(aFieldLinearTimeInterpolate, giveEvaluationOnTheInterpolateTimeForLinear) { - double interpolateTime = 0.2; + floater_t<4> interpolateTime = 0.2; this->destNew->setTime(interpolateTime); auto& layout = this->srcOld->gridLayout; diff --git a/tests/amr/data/field/variable/test_field_variable.cpp b/tests/amr/data/field/variable/test_field_variable.cpp index 5c6e52aaf..73be2fc0c 100644 --- a/tests/amr/data/field/variable/test_field_variable.cpp +++ b/tests/amr/data/field/variable/test_field_variable.cpp @@ -50,7 +50,7 @@ struct FieldVariableTest : public ::testing::TestWithParam using FieldVar = FieldVariable>, - Grid, HybridQuantity::Scalar>>; + Grid>, HybridQuantity::Scalar>>; // The interp order is of no importance to know if a quantity diff --git a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp index 86ce92e6d..c2106d8d7 100644 --- a/tests/amr/data/particles/copy/test_particledata_copyNd.cpp +++ b/tests/amr/data/particles/copy/test_particledata_copyNd.cpp @@ -115,12 +115,18 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy) this->sourceData.domainParticles.push_back(this->particle); this->destData.copy(this->sourceData); - EXPECT_THAT(this->destData.domainParticles[0].v, Pointwise(DoubleEq(), this->particle.v)); - EXPECT_THAT(this->destData.domainParticles[0].iCell, Eq(this->particle.iCell)); - EXPECT_THAT(this->destData.domainParticles[0].delta, - Pointwise(DoubleEq(), this->particle.delta)); - EXPECT_THAT(this->destData.domainParticles[0].weight, DoubleEq(this->particle.weight)); - EXPECT_THAT(this->destData.domainParticles[0].charge, DoubleEq(this->particle.charge)); + + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(this->destData.domainParticles[0].v, Pointwise(DoubleEq(), this->particle.v)); + EXPECT_THAT(this->destData.domainParticles[0].iCell, Eq(this->particle.iCell)); + EXPECT_THAT(this->destData.domainParticles[0].delta, + Pointwise(DoubleEq(), this->particle.delta)); + + + EXPECT_THAT(this->destData.domainParticles[0].weight, DoubleEq(this->particle.weight)); + EXPECT_THAT(this->destData.domainParticles[0].charge, DoubleEq(this->particle.charge)); + } // particle is in the domain of the source patchdata // and in last ghost of the destination patchdata @@ -130,12 +136,16 @@ TYPED_TEST(AParticlesDataND, PreservesAllParticleAttributesAfterCopy) this->sourceData.domainParticles.push_back(this->particle); this->destData.copy(this->sourceData); - EXPECT_THAT(this->destData.patchGhostParticles[0].v, Pointwise(DoubleEq(), this->particle.v)); - EXPECT_THAT(this->destData.patchGhostParticles[0].iCell, Eq(this->particle.iCell)); - EXPECT_THAT(this->destData.patchGhostParticles[0].delta, - Pointwise(DoubleEq(), this->particle.delta)); - EXPECT_THAT(this->destData.patchGhostParticles[0].weight, DoubleEq(this->particle.weight)); - EXPECT_THAT(this->destData.patchGhostParticles[0].charge, DoubleEq(this->particle.charge)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(this->destData.patchGhostParticles[0].v, + Pointwise(DoubleEq(), this->particle.v)); + EXPECT_THAT(this->destData.patchGhostParticles[0].iCell, Eq(this->particle.iCell)); + EXPECT_THAT(this->destData.patchGhostParticles[0].delta, + Pointwise(DoubleEq(), this->particle.delta)); + EXPECT_THAT(this->destData.patchGhostParticles[0].weight, DoubleEq(this->particle.weight)); + EXPECT_THAT(this->destData.patchGhostParticles[0].charge, DoubleEq(this->particle.charge)); + } } diff --git a/tests/amr/models/test_models.cpp b/tests/amr/models/test_models.cpp index 7c902de4e..d1c72eea7 100644 --- a/tests/amr/models/test_models.cpp +++ b/tests/amr/models/test_models.cpp @@ -21,6 +21,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include using namespace PHARE::core; @@ -33,13 +34,14 @@ using namespace PHARE::initializer::test_fn::func_1d; // density/etc are here static constexpr std::size_t dim = 1; static constexpr std::size_t interpOrder = 1; -using Field_t = Field<1, HybridQuantity::Scalar>; -using Grid1D = Grid, HybridQuantity::Scalar>; -using VecField1D = VecField; -using SymTensorField1D = SymTensorField; -using GridImplYee1D = GridLayoutImplYee; -using ParticleArray1D = ParticleArray; -using GridYee1D = GridLayout; +using core_Types = PHARE::core::PHARE_Types; + +using Field_t = core_Types::Field_t; +using Grid1D = core_Types::Grid_t; +using VecField1D = core_Types::VecField_t; +using SymTensorField1D = core_Types::SymTensorField_t; +using ParticleArray1D = core_Types::ParticleArray_t; +using GridYee1D = core_Types::GridLayout_t; using MaxwellianParticleInitializer1D = MaxwellianParticleInitializer; diff --git a/tests/amr/resources_manager/test_resources_manager.cpp b/tests/amr/resources_manager/test_resources_manager.cpp index 8ac7ff0d4..ce5011b36 100644 --- a/tests/amr/resources_manager/test_resources_manager.cpp +++ b/tests/amr/resources_manager/test_resources_manager.cpp @@ -27,20 +27,30 @@ using namespace PHARE::initializer::test_fn::func_1d; // density/etc are here static constexpr std::size_t dim = 1; static constexpr std::size_t interpOrder = 1; -using GridImplYee1D = GridLayoutImplYee; -using GridYee1D = GridLayout; - -using GridImplYee1D = GridLayoutImplYee; -using GridYee1D = GridLayout; -using Field1D = Field; -using Grid1D = Grid, HybridQuantity::Scalar>; -using VecField1D = VecField; -using SymTensorField1D = SymTensorField; -using IonPopulation1D = IonPopulation, VecField1D, SymTensorField1D>; -using Ions1D = Ions; -using Electromag1D = Electromag; -using Electrons1D = Electrons; -using HybridState1D = HybridState; + +using core_Types = PHARE::core::PHARE_Types; + +using Field1D = core_Types::Field_t; +using Grid1D = core_Types::Grid_t; +using VecField1D = core_Types::VecField_t; +using SymTensorField1D = core_Types::SymTensorField_t; +// using ParticleArray1D = core_Types::ParticleArray_t; +using GridYee1D = core_Types::GridLayout_t; + +// using GridImplYee1D = GridLayoutImplYee; +// using GridYee1D = GridLayout; + +// using GridImplYee1D = GridLayoutImplYee; +// using GridYee1D = GridLayout; +// using Field1D = Field; +// using Grid1D = Grid>, HybridQuantity::Scalar>; +// using VecField1D = VecField; +// using SymTensorField1D = SymTensorField; +using IonPopulation1D = IonPopulation, VecField1D, SymTensorField1D>; +using Ions1D = Ions; +using Electromag1D = Electromag; +using Electrons1D = Electrons; +using HybridState1D = HybridState; using MaxwellianParticleInitializer1D = MaxwellianParticleInitializer, GridYee1D>; diff --git a/tests/amr/resources_manager/test_resources_manager.hpp b/tests/amr/resources_manager/test_resources_manager.hpp index c741113e5..e1bd775e9 100644 --- a/tests/amr/resources_manager/test_resources_manager.hpp +++ b/tests/amr/resources_manager/test_resources_manager.hpp @@ -1,9 +1,10 @@ #ifndef PHARE_TESTS_AMR_TOOLS_RESSOURCE_RESSOURCE_TEST_1D_HPP #define PHARE_TESTS_AMR_TOOLS_RESSOURCE_RESSOURCE_TEST_1D_HPP -#include +#include "phare_core.hpp" + #include "test_resources_manager_basic_hierarchy.hpp" #include "core/data/grid/grid.hpp" #include "core/data/grid/gridlayout.hpp" @@ -19,6 +20,8 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include + using namespace PHARE::core; using namespace PHARE::amr; @@ -29,9 +32,12 @@ template class aResourceUserCollection : public ::testing::Test { public: - using Grid_t = Grid, HybridQuantity::Scalar>; + using core_Types = PHARE::core::PHARE_Types<1, 1>; + using Grid_t = core_Types::Grid_t; + using GridLayout_t = core_Types::GridLayout_t; + std::unique_ptr hierarchy; - ResourcesManager>, Grid_t> resourcesManager; + ResourcesManager resourcesManager; ResourcesUsers users; diff --git a/tests/amr/tagging/test_tagging.cpp b/tests/amr/tagging/test_tagging.cpp index 238126749..502d364ee 100644 --- a/tests/amr/tagging/test_tagging.cpp +++ b/tests/amr/tagging/test_tagging.cpp @@ -25,9 +25,9 @@ TEST(test_tagger, fromFactoryValid) { using phare_types = PHARE::PHARE_Types<1, 1, 2>; PHARE::initializer::PHAREDict dict; - dict["model"] = std::string{"HybridModel"}; - dict["method"] = std::string{"default"}; - dict["threshold"] = 0.2; + dict["model"] = std::string{"HybridModel"}; + dict["method"] = std::string{"default"}; + dict["threshold"] = 0.2; auto hybridTagger = TaggerFactory::make(dict); EXPECT_TRUE(hybridTagger != nullptr); } @@ -36,23 +36,23 @@ TEST(test_tagger, fromFactoryInvalid) { using phare_types = PHARE::PHARE_Types<1, 1, 2>; PHARE::initializer::PHAREDict dict; - dict["model"] = std::string{"invalidModel"}; - dict["method"] = std::string{"invalidStrat"}; + dict["model"] = std::string{"invalidModel"}; + dict["method"] = std::string{"invalidStrat"}; auto hybridTagger = TaggerFactory::make(dict); - auto badTagger = TaggerFactory::make(dict); + auto badTagger = TaggerFactory::make(dict); EXPECT_TRUE(badTagger == nullptr); } -using Param = std::vector; -using RetType = std::shared_ptr>; +using Param = std::vector>; +using RetType = std::shared_ptr>>; RetType step1(Param const& x) { - std::vector values(x.size()); + std::vector> values(x.size()); std::transform(std::begin(x), std::end(x), std::begin(values), - [](auto xx) { return std::tanh((xx - 0.52) / 0.05); }); - return std::make_shared>(std::move(values)); + [](auto xx) { return std::tanh((xx - 0.52f) / 0.05f); }); + return std::make_shared>>(std::move(values)); } RetType step2(Param const& x, Param const& y) diff --git a/tests/core/data/electrons/test_electrons.cpp b/tests/core/data/electrons/test_electrons.cpp index d69df65ff..368853dd7 100644 --- a/tests/core/data/electrons/test_electrons.cpp +++ b/tests/core/data/electrons/test_electrons.cpp @@ -133,8 +133,8 @@ struct ElectronsTest : public ::testing::Test using GridYee = GridLayout>; - using GridND = Grid, HybridQuantity::Scalar>; - using FieldND = Field; + using GridND = Grid>, HybridQuantity::Scalar>; + using FieldND = Field>; using VecFieldND = VecField; using SymTensorFieldND = SymTensorField; using ParticleArray_t = ParticleArray; @@ -372,7 +372,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsDensityEqualIonDensity) for (std::uint32_t i = psi_X; i < pei_X; ++i) { - EXPECT_DOUBLE_EQ(Ni(i), Ne(i)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Ni(i), Ne(i)); + } } } else if constexpr (dim == 2) @@ -386,7 +389,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsDensityEqualIonDensity) { for (std::uint32_t j = psi_Y; j < pei_Y; ++j) { - EXPECT_DOUBLE_EQ(Ni(i, j), Ne(i, j)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Ni(i, j), Ne(i, j)); + } } } } @@ -405,7 +411,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsDensityEqualIonDensity) { for (std::uint32_t k = psi_Z; k < pei_Z; ++k) { - EXPECT_DOUBLE_EQ(Ni(i, j, k), Ne(i, j, k)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Ni(i, j, k), Ne(i, j, k)); + } } } } @@ -437,7 +446,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsVelocityEqualIonVelocityMinusJ) { auto const JOnV = GridYee::project(Jcomp, {i}, projector()); - EXPECT_DOUBLE_EQ(Vecomp(i), Vicomp(i) - JOnV / Ne_(i)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Vecomp(i), Vicomp(i) - JOnV / Ne_(i)); + } } } else if constexpr (dim == 2) @@ -453,7 +465,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsVelocityEqualIonVelocityMinusJ) { auto const JOnV = GridYee::project(Jcomp, {i, j}, projector()); - EXPECT_DOUBLE_EQ(Vecomp(i, j), Vicomp(i, j) - JOnV / Ne_(i, j)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Vecomp(i, j), Vicomp(i, j) - JOnV / Ne_(i, j)); + } } } } @@ -474,7 +489,11 @@ TYPED_TEST(ElectronsTest, ThatElectronsVelocityEqualIonVelocityMinusJ) { auto const JOnV = GridYee::project(Jcomp, {i, j, k}, projector()); - EXPECT_DOUBLE_EQ(Vecomp(i, j, k), Vicomp(i, j, k) - JOnV / Ne_(i, j, k)); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Vecomp(i, j, k), + Vicomp(i, j, k) - JOnV / Ne_(i, j, k)); + } } } } @@ -510,7 +529,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsPressureEqualsNeTe) for (std::uint32_t i = psi_X; i < pei_X; ++i) { - EXPECT_DOUBLE_EQ(Pe_(i), Ne_(i) * Te); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Pe_(i), Ne_(i) * Te); + } } } else if constexpr (dim == 2) @@ -524,7 +546,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsPressureEqualsNeTe) { for (std::uint32_t j = psi_Y; j < pei_Y; ++j) { - EXPECT_DOUBLE_EQ(Pe_(i, j), Ne_(i, j) * Te); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Pe_(i, j), Ne_(i, j) * Te); + } } } } @@ -543,7 +568,10 @@ TYPED_TEST(ElectronsTest, ThatElectronsPressureEqualsNeTe) { for (std::uint32_t k = psi_Z; k < pei_Z; ++k) { - EXPECT_DOUBLE_EQ(Pe_(i, j, k), Ne_(i, j, k) * Te); + if constexpr (std::is_same_v, double>) + { + EXPECT_DOUBLE_EQ(Pe_(i, j, k), Ne_(i, j, k) * Te); + } } } } diff --git a/tests/core/data/field/test_field_fixtures.hpp b/tests/core/data/field/test_field_fixtures.hpp index c63c1ce97..24ca43034 100644 --- a/tests/core/data/field/test_field_fixtures.hpp +++ b/tests/core/data/field/test_field_fixtures.hpp @@ -1,13 +1,14 @@ #ifndef PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP #define PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP +#include "core/def.hpp" #include "core/data/field/field.hpp" namespace PHARE::core { template -using Field_t = Field; +using Field_t = Field>; } // namespace PHARE::core diff --git a/tests/core/data/gridlayout/gridlayout_allocsize.hpp b/tests/core/data/gridlayout/gridlayout_allocsize.hpp index 0471f6da2..1d27e546d 100644 --- a/tests/core/data/gridlayout/gridlayout_allocsize.hpp +++ b/tests/core/data/gridlayout/gridlayout_allocsize.hpp @@ -76,14 +76,14 @@ auto createAllocSizeParam() while (inputFile >> iQuantity) { std::array numberCells; - std::array dl; + std::array, GridLayoutImpl::dimension> dl; writeToArray(inputFile, numberCells); writeToArray(inputFile, dl); params.emplace_back(); params.back().base = createParam( - dl, numberCells, Point{}); + dl, numberCells, Point, GridLayoutImpl::dimension>{}); writeToArray(inputFile, params.back().expectedAllocSize); writeToArray(inputFile, params.back().expectedAllocSizeDerived); diff --git a/tests/core/data/gridlayout/gridlayout_base_params.hpp b/tests/core/data/gridlayout/gridlayout_base_params.hpp index c34752134..ddde1e0cc 100644 --- a/tests/core/data/gridlayout/gridlayout_base_params.hpp +++ b/tests/core/data/gridlayout/gridlayout_base_params.hpp @@ -22,10 +22,10 @@ struct GridLayoutTestParam using Grid_t = Grid{})), HybridQuantity::Scalar>; std::shared_ptr> layout; - std::array dxdydz; + std::array, dim> dxdydz; std::array nbCellXYZ; - Point origin; + Point, dim> origin; HybridQuantity::Scalar currentQuantity; diff --git a/tests/core/data/gridlayout/gridlayout_cell_centered_coord.hpp b/tests/core/data/gridlayout/gridlayout_cell_centered_coord.hpp index 652f39c19..9112c494b 100644 --- a/tests/core/data/gridlayout/gridlayout_cell_centered_coord.hpp +++ b/tests/core/data/gridlayout/gridlayout_cell_centered_coord.hpp @@ -19,9 +19,9 @@ struct GridLayoutCellCenteringParam GridLayoutTestParam base; std::vector> iCellForCentering; - std::vector> expectedPosition; + std::vector, GridLayoutImpl::dimension>> expectedPosition; - std::vector> actualPosition; + std::vector, GridLayoutImpl::dimension>> actualPosition; template auto cellCenteredCoord_impl(Array const& array, std::index_sequence) @@ -45,10 +45,10 @@ struct GridLayoutCellCenteringParam for (auto&& iCell : iCellForCentering) { - Point pos; + Point, GridLayoutImpl::dimension> pos; pos = cellCenteredCoord(iCell); - std::array actualPos; + std::array, GridLayoutImpl::dimension> actualPos; for (std::size_t iDim = 0; iDim < GridLayoutImpl::dimension; ++iDim) { @@ -92,12 +92,12 @@ auto createCellCenteringParam() {"Vz", HybridQuantity::Scalar::Vz}, {"P", HybridQuantity::Scalar::P}}; std::array nbCell; - std::array dl; + std::array, GridLayoutImpl::dimension> dl; std::array iStart; std::array iEnd; - std::array origin; + std::array, GridLayoutImpl::dimension> origin; writeToArray(summary, nbCell); writeToArray(summary, dl); @@ -110,12 +110,12 @@ auto createCellCenteringParam() params.emplace_back(); // NOTE: c++17 : Point{origin}, C++14 : Point{origin} - params.back().base - = createParam(dl, nbCell, Point{origin}); + params.back().base = createParam( + dl, nbCell, Point, GridLayoutImpl::dimension>{origin}); std::array icell; - std::array realPosition; + std::array, GridLayoutImpl::dimension> realPosition; while (writeToArray(value, icell) && writeToArray(value, realPosition)) { diff --git a/tests/core/data/gridlayout/gridlayout_deriv.cpp b/tests/core/data/gridlayout/gridlayout_deriv.cpp index 4a598e5ef..94af6f5e6 100644 --- a/tests/core/data/gridlayout/gridlayout_deriv.cpp +++ b/tests/core/data/gridlayout/gridlayout_deriv.cpp @@ -1,12 +1,15 @@ -#include +// #include #include "gridlayout_deriv.hpp" #include "gridlayout_params.hpp" #include "gridlayout_test.hpp" +#include -std::vector read(std::string filename) +PHARE::core::floater_t<4> constexpr static _PI_ = std::numbers::pi_v>; + +auto read(std::string filename) { std::ifstream readFile(filename); assert(readFile.is_open()); @@ -14,7 +17,18 @@ std::vector read(std::string filename) std::copy(std::istream_iterator(readFile), std::istream_iterator(), std::back_inserter(x)); - return x; + + if constexpr (std::is_same_v, double>) + { + return x; + } + else + { + std::vector> fx(x.size()); + for (std::size_t i = 0; i < x.size(); i++) + fx[i] = x[i]; + return fx; + } } // ----------------------------------------------------------------------------- @@ -42,14 +56,17 @@ TYPED_TEST(a1DDerivative, DXBY1D) for (auto ix = 0u; ix <= gei_d_X; ++ix) { - Point point = this->layout.fieldNodeCoordinates(this->By, {0.}, ix); - this->By(ix) = std::cos(2 * M_PI / 5. * point[0]); + auto const point = this->layout.fieldNodeCoordinates(this->By, {0.}, ix); + this->By(ix) = std::cos(2 * _PI_ / 5.f * point[0]); } for (auto ix = psi_p_X; ix <= pei_p_X; ++ix) { auto localDerivative = this->layout.template deriv(this->By, {ix}); - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[ix], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[ix], 1e-12)); + } } } @@ -68,14 +85,17 @@ TYPED_TEST(a1DDerivative, DXEZ1D) for (auto ix = 0u; ix <= gei_p_X; ++ix) { - Point point = this->layout.fieldNodeCoordinates(this->Ez, {0.}, ix); - this->Ez(ix) = std::cos(2 * M_PI / 5. * point[0]); + auto const point = this->layout.fieldNodeCoordinates(this->Ez, {0.}, ix); + this->Ez(ix) = std::cos(2 * _PI_ / 5.f * point[0]); } for (auto ix = psi_d_X; ix <= pei_d_X; ++ix) { auto localDerivative = this->layout.template deriv(this->Ez, {ix}); - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[ix], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[ix], 1e-12)); + } } } @@ -109,9 +129,9 @@ TYPED_TEST(a2DDerivative, DXBY2D) { for (auto iy = 0u; iy <= gei_p_Y; ++iy) { - Point point = this->layout.fieldNodeCoordinates(this->By, {0., 0.}, ix, iy); + auto const point = this->layout.fieldNodeCoordinates(this->By, {0., 0.}, ix, iy); this->By(ix, iy) - = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -124,7 +144,10 @@ TYPED_TEST(a2DDerivative, DXBY2D) { auto localDerivative = this->layout.template deriv(this->By, {ix, iy}); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -148,9 +171,9 @@ TYPED_TEST(a2DDerivative, DYBY2D) { for (auto iy = 0u; iy <= gei_p_Y; ++iy) { - Point point = this->layout.fieldNodeCoordinates(this->By, {0., 0.}, ix, iy); + auto const point = this->layout.fieldNodeCoordinates(this->By, {0., 0.}, ix, iy); this->By(ix, iy) - = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -162,7 +185,10 @@ TYPED_TEST(a2DDerivative, DYBY2D) { auto localDerivative = this->layout.template deriv(this->By, {ix, iy}); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -186,9 +212,9 @@ TYPED_TEST(a2DDerivative, DXEZ2D) { for (auto iy = 0u; iy <= gei_p_Y; ++iy) { - Point point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0.}, ix, iy); + auto const point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0.}, ix, iy); this->Ez(ix, iy) - = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -200,7 +226,10 @@ TYPED_TEST(a2DDerivative, DXEZ2D) { auto localDerivative = this->layout.template deriv(this->Ez, {ix, iy}); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -225,7 +254,7 @@ TYPED_TEST(a2DDerivative, DYEZ2D) { auto point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0.}, ix, iy); this->Ez(ix, iy) - = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -237,7 +266,10 @@ TYPED_TEST(a2DDerivative, DYEZ2D) { auto localDerivative = this->layout.template deriv(this->Ez, {ix, iy}); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -278,9 +310,9 @@ TYPED_TEST(a3DDerivative, DXBY3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->By, {0., 0., 0.}, ix, iy, iz); - this->By(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->By(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -296,7 +328,10 @@ TYPED_TEST(a3DDerivative, DXBY3D) auto localDerivative = this->layout.template deriv(this->By, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -327,9 +362,9 @@ TYPED_TEST(a3DDerivative, DYBY3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->By, {0., 0., 0.}, ix, iy, iz); - this->By(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->By(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -345,7 +380,10 @@ TYPED_TEST(a3DDerivative, DYBY3D) auto localDerivative = this->layout.template deriv(this->By, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -377,9 +415,9 @@ TYPED_TEST(a3DDerivative, DZBY3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->By, {0., 0., 0.}, ix, iy, iz); - this->By(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->By(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -395,7 +433,10 @@ TYPED_TEST(a3DDerivative, DZBY3D) auto localDerivative = this->layout.template deriv(this->By, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -426,9 +467,9 @@ TYPED_TEST(a3DDerivative, DXEZ3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0., 0.}, ix, iy, iz); - this->Ez(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->Ez(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -444,7 +485,10 @@ TYPED_TEST(a3DDerivative, DXEZ3D) auto localDerivative = this->layout.template deriv(this->Ez, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -475,9 +519,9 @@ TYPED_TEST(a3DDerivative, DYEZ3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0., 0.}, ix, iy, iz); - this->Ez(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->Ez(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -493,7 +537,10 @@ TYPED_TEST(a3DDerivative, DYEZ3D) auto localDerivative = this->layout.template deriv(this->Ez, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } @@ -524,9 +571,9 @@ TYPED_TEST(a3DDerivative, DZEZ3D) for (auto iz = 0u; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates(this->Ez, {0., 0., 0.}, ix, iy, iz); - this->Ez(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + this->Ez(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -542,7 +589,10 @@ TYPED_TEST(a3DDerivative, DZEZ3D) auto localDerivative = this->layout.template deriv(this->Ez, {ix, iy, iz}); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localDerivative, ::testing::DoubleNear(expDerValue[index_], 1e-12)); + } } } } diff --git a/tests/core/data/gridlayout/gridlayout_deriv.hpp b/tests/core/data/gridlayout/gridlayout_deriv.hpp index 9cf6e4b2e..6fc631a04 100644 --- a/tests/core/data/gridlayout/gridlayout_deriv.hpp +++ b/tests/core/data/gridlayout/gridlayout_deriv.hpp @@ -17,7 +17,7 @@ using namespace PHARE::core; -std::vector read(std::string filename); +auto read(std::string filename); diff --git a/tests/core/data/gridlayout/gridlayout_field_centered_coord.hpp b/tests/core/data/gridlayout/gridlayout_field_centered_coord.hpp index e7a363012..6f5e05eb3 100644 --- a/tests/core/data/gridlayout/gridlayout_field_centered_coord.hpp +++ b/tests/core/data/gridlayout/gridlayout_field_centered_coord.hpp @@ -17,8 +17,8 @@ struct GridLayoutFieldCenteringParam GridLayoutTestParam base; std::vector> iCellForCentering; - std::vector> expectedPosition; - std::vector> actualPosition; + std::vector, GridLayoutImpl::dimension>> expectedPosition; + std::vector, GridLayoutImpl::dimension>> actualPosition; template auto fieldCoord_impl(Array const& array, std::index_sequence) @@ -49,7 +49,7 @@ struct GridLayoutFieldCenteringParam Point pos; pos = fieldCoord(iCell); - std::array actualPos; + std::array, GridLayoutImpl::dimension> actualPos; for (std::size_t iDim = 0; iDim < GridLayoutImpl::dimension; ++iDim) { @@ -99,10 +99,10 @@ auto createFieldCenteringParam() while (summary >> quantity) { std::array nbCell; - std::array dl; + std::array, GridLayoutImpl::dimension> dl; std::array iStart; std::array iEnd; - std::array origin; + std::array, GridLayoutImpl::dimension> origin; writeToArray(summary, nbCell); writeToArray(summary, dl); @@ -123,7 +123,7 @@ auto createFieldCenteringParam() std::array icell; - std::array realPosition; + std::array, GridLayoutImpl::dimension> realPosition; while (value >> quantity && writeToArray(value, icell) && writeToArray(value, realPosition)) { diff --git a/tests/core/data/gridlayout/gridlayout_indexing.hpp b/tests/core/data/gridlayout/gridlayout_indexing.hpp index 6f3223f25..23962f490 100644 --- a/tests/core/data/gridlayout/gridlayout_indexing.hpp +++ b/tests/core/data/gridlayout/gridlayout_indexing.hpp @@ -85,14 +85,14 @@ auto createIndexingParam() while (inputFile >> iQuantity) { std::array numberCells; - std::array dl; + std::array, GridLayoutImpl::dimension> dl; writeToArray(inputFile, numberCells); writeToArray(inputFile, dl); params.emplace_back(); params.back().base = createParam( - dl, numberCells, Point{}); + dl, numberCells, Point, GridLayoutImpl::dimension>{}); writeToArray(inputFile, params.back().expectedPSI); writeToArray(inputFile, params.back().expectedPEI); diff --git a/tests/core/data/gridlayout/gridlayout_laplacian.cpp b/tests/core/data/gridlayout/gridlayout_laplacian.cpp index 24ee8555c..0c7afe948 100644 --- a/tests/core/data/gridlayout/gridlayout_laplacian.cpp +++ b/tests/core/data/gridlayout/gridlayout_laplacian.cpp @@ -32,13 +32,20 @@ TYPED_TEST(a1DLaplacian, LaplacianJx1D) for (auto ix = gsi_X; ix <= gei_X; ++ix) { auto point = this->layout.fieldNodeCoordinates(this->Jx, Point{0.}, ix); - this->Jx(ix) = std::sinh(0.1 * point[0]); + this->Jx(ix) = std::sinh(0.1f * point[0]); } for (auto ix = psi_X; ix <= pei_X; ++ix) { auto localLaplacian = this->layout.laplacian(this->Jx, make_index(ix)); - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + } + else + { + EXPECT_THAT(localLaplacian, ::testing::FloatNear(expLapValue[ix], 1e-6)); + } } } @@ -60,13 +67,16 @@ TYPED_TEST(a1DLaplacian, LaplacianJy1D) for (auto ix = gsi_X; ix <= gei_X; ++ix) { auto point = this->layout.fieldNodeCoordinates(this->Jy, Point{0.}, ix); - this->Jy(ix) = std::sinh(0.3 * point[0]); + this->Jy(ix) = std::sinh(0.3f * point[0]); } for (auto ix = psi_X; ix <= pei_X; ++ix) { auto localLaplacian = this->layout.laplacian(this->Jy, make_index(ix)); - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + } } } @@ -88,13 +98,16 @@ TYPED_TEST(a1DLaplacian, LaplacianJz1D) for (auto ix = gsi_X; ix <= gei_X; ++ix) { auto point = this->layout.fieldNodeCoordinates(this->Jz, Point{0.}, ix); - this->Jz(ix) = std::sinh(0.2 * point[0]); + this->Jz(ix) = std::sinh(0.2f * point[0]); } for (std::uint32_t ix = psi_X; ix <= pei_X; ++ix) { auto localLaplacian = this->layout.laplacian(this->Jz, make_index(ix)); - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[ix], 1e-12)); + } } } @@ -132,7 +145,7 @@ TYPED_TEST(a2DLaplacian, LaplacianJx2D) for (auto iy = gsi_Y; iy <= gei_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(this->Jx, Point{0., 0.}, ix, iy); - this->Jx(ix, iy) = std::sinh(0.1 * point[0]) * std::cosh(0.1 * point[1]); + this->Jx(ix, iy) = std::sinh(0.1f * point[0]) * std::cosh(0.1f * point[1]); } } @@ -146,7 +159,10 @@ TYPED_TEST(a2DLaplacian, LaplacianJx2D) auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } @@ -174,7 +190,7 @@ TYPED_TEST(a2DLaplacian, LaplacianJy2D) for (auto iy = gsi_Y; iy <= gei_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(this->Jy, Point{0., 0.}, ix, iy); - this->Jy(ix, iy) = std::sinh(0.3 * point[0]) * std::cosh(0.3 * point[1]); + this->Jy(ix, iy) = std::sinh(0.3f * point[0]) * std::cosh(0.3f * point[1]); } } @@ -188,7 +204,10 @@ TYPED_TEST(a2DLaplacian, LaplacianJy2D) auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } @@ -216,7 +235,7 @@ TYPED_TEST(a2DLaplacian, LaplacianJz2D) for (auto iy = gsi_Y; iy <= gei_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(this->Jz, Point{0., 0.}, ix, iy); - this->Jz(ix, iy) = std::sinh(0.2 * point[0]) * std::cosh(0.2 * point[1]); + this->Jz(ix, iy) = std::sinh(0.2f * point[0]) * std::cosh(0.2f * point[1]); } } @@ -230,7 +249,10 @@ TYPED_TEST(a2DLaplacian, LaplacianJz2D) auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } @@ -277,8 +299,8 @@ TYPED_TEST(a3DLaplacian, LaplacianJx3D) { auto point = this->layout.fieldNodeCoordinates(this->Jx, Point{0., 0., 0.}, ix, iy, iz); - this->Jx(ix, iy, iz) = std::sinh(0.1 * point[0]) * std::cosh(0.1 * point[1]) - * std::tanh(0.1 * point[2]); + this->Jx(ix, iy, iz) = std::sinh(0.1f * point[0]) * std::cosh(0.1f * point[1]) + * std::tanh(0.1f * point[2]); } } } @@ -295,7 +317,10 @@ TYPED_TEST(a3DLaplacian, LaplacianJx3D) auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } @@ -332,8 +357,8 @@ TYPED_TEST(a3DLaplacian, LaplacianJy3D) { auto point = this->layout.fieldNodeCoordinates(this->Jy, Point{0., 0., 0.}, ix, iy, iz); - this->Jy(ix, iy, iz) = std::sinh(0.3 * point[0]) * std::cosh(0.3 * point[1]) - * std::tanh(0.3 * point[2]); + this->Jy(ix, iy, iz) = std::sinh(0.3f * point[0]) * std::cosh(0.3f * point[1]) + * std::tanh(0.3f * point[2]); } } } @@ -350,7 +375,10 @@ TYPED_TEST(a3DLaplacian, LaplacianJy3D) auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } @@ -387,8 +415,8 @@ TYPED_TEST(a3DLaplacian, LaplacianJz3D) { auto point = this->layout.fieldNodeCoordinates(this->Jz, Point{0., 0., 0.}, ix, iy, iz); - this->Jz(ix, iy, iz) = std::sinh(0.2 * point[0]) * std::cosh(0.2 * point[1]) - * std::tanh(0.2 * point[2]); + this->Jz(ix, iy, iz) = std::sinh(0.2f * point[0]) * std::cosh(0.2f * point[1]) + * std::tanh(0.2f * point[2]); } } } @@ -405,7 +433,10 @@ TYPED_TEST(a3DLaplacian, LaplacianJz3D) auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(localLaplacian, ::testing::DoubleNear(expLapValue[index_], 1e-12)); + } } } } diff --git a/tests/core/data/gridlayout/gridlayout_params.hpp b/tests/core/data/gridlayout/gridlayout_params.hpp index 92f0aa2ec..aaa15bc8a 100644 --- a/tests/core/data/gridlayout/gridlayout_params.hpp +++ b/tests/core/data/gridlayout/gridlayout_params.hpp @@ -51,9 +51,9 @@ inline HybridQuantity::Scalar getQuantity(std::uint32_t iQuantity) template -auto createParam(std::array const& dxdydz, +auto createParam(std::array, GridLayoutImpl::dimension> const& dxdydz, std::array const& nbCellXYZ, - Point const& origin) + Point, GridLayoutImpl::dimension> const& origin) { GridLayoutTestParam param{}; param.layout = std::make_shared>(dxdydz, nbCellXYZ, origin); diff --git a/tests/core/data/gridlayout/test_gridlayout.hpp b/tests/core/data/gridlayout/test_gridlayout.hpp index d31f630a6..5d1271786 100644 --- a/tests/core/data/gridlayout/test_gridlayout.hpp +++ b/tests/core/data/gridlayout/test_gridlayout.hpp @@ -15,9 +15,10 @@ class TestGridLayout : public GridLayout TestGridLayout() = default; TestGridLayout(std::uint32_t cells) - : GridLayout{PHARE::core::ConstArray(1.0 / cells), + : GridLayout{PHARE::core::ConstArray, dim>(1.0 / cells), PHARE::core::ConstArray(cells), - PHARE::core::Point{PHARE::core::ConstArray(0)}} + PHARE::core::Point, dim>{ + PHARE::core::ConstArray, dim>(0)}} { } diff --git a/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp b/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp index f3f4899cb..4fcbef90c 100644 --- a/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp +++ b/tests/core/data/maxwellian_particle_initializer/test_maxwellian_particle_initializer.cpp @@ -60,6 +60,8 @@ TEST_F(AMaxwellianParticleInitializer1D, loadsTheCorrectNbrOfParticles) TEST_F(AMaxwellianParticleInitializer1D, loadsParticlesInTheDomain) { + floater_t<4> constexpr static ZRO = 0; + initializer->loadParticles(particles, layout); for (auto const& particle : particles) { @@ -67,9 +69,9 @@ TEST_F(AMaxwellianParticleInitializer1D, loadsParticlesInTheDomain) auto pos = positionAsPoint(particle, layout); auto endDomain = layout.origin()[0] + layout.nbrCells()[0] * layout.meshSize()[0]; - if (!((pos[0] > 0.) and (pos[0] < endDomain))) + if (!((pos[0] > ZRO) and (pos[0] < endDomain))) std::cout << "position : " << pos[0] << " not in domain (0," << endDomain << ")\n"; - EXPECT_TRUE(pos[0] > 0. && pos[0] < endDomain); + EXPECT_TRUE(pos[0] > ZRO && pos[0] < endDomain); } } diff --git a/tests/core/data/particles/test_interop.cpp b/tests/core/data/particles/test_interop.cpp index c9eb98f12..4714e99f1 100644 --- a/tests/core/data/particles/test_interop.cpp +++ b/tests/core/data/particles/test_interop.cpp @@ -32,8 +32,8 @@ TYPED_TEST(ParticleListTest, SoAandAoSInterop) view.weight = 1 + i; view.charge = 1 + i; view.iCell = ConstArray(i); - view.delta = ConstArray(i + 1); - view.v = ConstArray(view.weight + 2); + view.delta = ConstArray, dim>(i + 1); + view.v = ConstArray, 3>(view.weight + 2); EXPECT_EQ(std::copy(view), view); } EXPECT_EQ(contiguous.size(), size); diff --git a/tests/core/data/particles/test_main.cpp b/tests/core/data/particles/test_main.cpp index f3ad991c1..b9214f93a 100644 --- a/tests/core/data/particles/test_main.cpp +++ b/tests/core/data/particles/test_main.cpp @@ -74,8 +74,8 @@ TEST_F(AParticle, CanBeReducedToAnIntegerPoint) TEST_F(AParticle, CanBeReducedToAnAbsolutePositionPoint) { - Point origin; - std::array meshSize{{0.2, 0.05, 0.4}}; + Point, 3> origin; + std::array, 3> meshSize{{0.2, 0.05, 0.4}}; std::array nbrCells{{20, 30, 40}}; GridLayout> layout{meshSize, nbrCells, origin, Box{Point{40, 60, 80}, Point{59, 89, 119}}}; @@ -83,12 +83,12 @@ TEST_F(AParticle, CanBeReducedToAnAbsolutePositionPoint) auto iCell = layout.AMRToLocal(Point{part.iCell}); auto p = positionAsPoint(part, layout); auto startIndexes = layout.physicalStartIndex(QtyCentering::primal); - auto expectedPosition = Point{}; + auto expectedPosition = Point, 3>{}; for (auto i = 0u; i < 3; ++i) { expectedPosition[i] - = origin[i] + meshSize[i] * (iCell[i] - startIndexes[i] + part.delta[i]); - EXPECT_DOUBLE_EQ(expectedPosition[i], p[i]); + = origin[i] + float(meshSize[i]) * (iCell[i] - startIndexes[i] + part.delta[i]); + EXPECT_FLOAT_EQ(expectedPosition[i], p[i]); } } diff --git a/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp b/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp index cb20297dc..85700a2e7 100644 --- a/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp +++ b/tests/core/data/tensorfield/test_tensorfield_fixtures.hpp @@ -23,8 +23,8 @@ class UsableTensorField : public TensorField, HybridQuantity, rank_ public: auto static constexpr dimension = dim; using Super = TensorField, HybridQuantity, rank_>; - using Grid_t = Grid, HybridQuantity::Scalar>; - using tensor_t = typename Super::tensor_t; + using Grid_t = Grid>, HybridQuantity::Scalar>; + using tensor_t = typename Super::tensor_t; template UsableTensorField(std::string const& name, GridLayout const& layout, tensor_t qty) @@ -50,9 +50,8 @@ class UsableTensorField : public TensorField, HybridQuantity, rank_ auto static make_grids(ComponentNames const& compNames, GridLayout const& layout, tensor_t qty) { auto qts = HybridQuantity::componentsQuantities(qty); - return for_N([&](auto i) { - return Grid_t{compNames[i], qts[i], layout.allocSize(qts[i])}; - }); + return for_N( + [&](auto i) { return Grid_t{compNames[i], qts[i], layout.allocSize(qts[i])}; }); } std::array xyz; diff --git a/tests/core/numerics/ampere/test_main.cpp b/tests/core/numerics/ampere/test_main.cpp index ec9827539..c8b546aa0 100644 --- a/tests/core/numerics/ampere/test_main.cpp +++ b/tests/core/numerics/ampere/test_main.cpp @@ -1,8 +1,6 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" -#include -#include #include "core/data/grid/grid.hpp" @@ -21,6 +19,13 @@ #include "tests/core/data/gridlayout/gridlayout_test.hpp" +#include +#include +#include + +PHARE::core::floater_t<4> constexpr static _PI_ = std::numbers::pi_v>; + + using namespace PHARE::core; @@ -206,8 +211,8 @@ TEST_F(Ampere1DTest, ampere1DCalculatedOk) for (std::uint32_t ix = gsi_d_X; ix <= gei_d_X; ++ix) { auto point = this->layout.fieldNodeCoordinates(By, Point{0.}, ix); - By(ix) = std::cos(2 * M_PI / 5. * point[0]); - Bz(ix) = std::sin(2 * M_PI / 5. * point[0]); + By(ix) = std::cos(2 * _PI_ / 5.f * point[0]); + Bz(ix) = std::sin(2 * _PI_ / 5.f * point[0]); } ampere.setLayout(&layout); @@ -218,8 +223,16 @@ TEST_F(Ampere1DTest, ampere1DCalculatedOk) for (std::uint32_t ix = psi_p_X; ix <= pei_p_X; ++ix) { - EXPECT_THAT(Jy(ix), ::testing::DoubleNear((expectedJy[ix]), 1e-12)); - EXPECT_THAT(Jz(ix), ::testing::DoubleNear((expectedJz[ix]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Jy(ix), ::testing::DoubleNear((expectedJy[ix]), 1e-12)); + EXPECT_THAT(Jz(ix), ::testing::DoubleNear((expectedJz[ix]), 1e-12)); + } + else + { + EXPECT_THAT(Jy(ix), ::testing::FloatNear((expectedJy[ix]), 1e-5)); + EXPECT_THAT(Jz(ix), ::testing::FloatNear((expectedJz[ix]), 1e-5)); + } } } @@ -250,7 +263,7 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (std::uint32_t iy = gsi_d_Y; iy <= gei_d_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(Bx, Point{0., 0.}, ix, iy); - Bx(ix, iy) = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + Bx(ix, iy) = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -259,7 +272,7 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (std::uint32_t iy = gsi_p_Y; iy <= gei_p_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(By, Point{0., 0.}, ix, iy); - By(ix, iy) = std::cos(2 * M_PI / 5. * point[0]) * std::tanh(2 * M_PI / 6. * point[1]); + By(ix, iy) = std::cos(2 * _PI_ / 5.f * point[0]) * std::tanh(2 * _PI_ / 6.f * point[1]); } } @@ -268,7 +281,7 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (std::uint32_t iy = gsi_d_Y; iy <= gei_d_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(Bz, Point{0., 0.}, ix, iy); - Bz(ix, iy) = std::sin(2 * M_PI / 5. * point[0]) * std::tanh(2 * M_PI / 6. * point[1]); + Bz(ix, iy) = std::sin(2 * _PI_ / 5.f * point[0]) * std::tanh(2 * _PI_ / 6.f * point[1]); } } @@ -291,7 +304,14 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (auto iy = psi_p_Y; iy <= pei_p_Y; ++iy) { std::uint32_t index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Jx(ix, iy), ::testing::DoubleNear((expectedJx[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Jx(ix, iy), ::testing::DoubleNear((expectedJx[index_]), 1e-12)); + } + else + { + EXPECT_THAT(Jx(ix, iy), ::testing::FloatNear((expectedJx[index_]), 1e-6)); + } } } @@ -302,7 +322,14 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (auto iy = psi_d_Y; iy <= pei_d_Y; ++iy) { std::uint32_t index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Jy(ix, iy), ::testing::DoubleNear((expectedJy[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Jy(ix, iy), ::testing::DoubleNear((expectedJy[index_]), 1e-12)); + } + else + { + EXPECT_THAT(Jy(ix, iy), ::testing::FloatNear((expectedJy[index_]), 1e-6)); + } } } @@ -313,7 +340,14 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) for (auto iy = psi_p_Y; iy <= pei_p_Y; ++iy) { std::uint32_t index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Jz(ix, iy), ::testing::DoubleNear((expectedJz[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Jz(ix, iy), ::testing::DoubleNear((expectedJz[index_]), 1e-12)); + } + else + { + EXPECT_THAT(Jz(ix, iy), ::testing::FloatNear((expectedJz[index_]), 1e-6)); + } } } } @@ -351,11 +385,11 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) { for (std::uint32_t iz = gsi_d_Z; iz <= gei_d_Z; ++iz) { - Point point = this->layout.fieldNodeCoordinates( + auto const point = this->layout.fieldNodeCoordinates( Bx, Point{0., 0., 0.}, ix, iy, iz); - Bx(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::tanh(2 * M_PI / 12. * point[2]); + Bx(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::tanh(2 * _PI_ / 12.f * point[2]); } } } @@ -366,11 +400,11 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) { for (std::uint32_t iz = gsi_d_Z; iz <= gei_d_Z; ++iz) { - Point point = this->layout.fieldNodeCoordinates( + auto const point = this->layout.fieldNodeCoordinates( By, Point{0., 0., 0.}, ix, iy, iz); - By(ix, iy, iz) = std::tanh(2 * M_PI / 5. * point[0]) - * std::sin(2 * M_PI / 6. * point[1]) - * std::cos(2 * M_PI / 12. * point[2]); + By(ix, iy, iz) = std::tanh(2 * _PI_ / 5.f * point[0]) + * std::sin(2 * _PI_ / 6.f * point[1]) + * std::cos(2 * _PI_ / 12.f * point[2]); } } } @@ -381,11 +415,11 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) { for (std::uint32_t iz = gsi_p_Z; iz <= gei_p_Z; ++iz) { - Point point = this->layout.fieldNodeCoordinates( + auto const point = this->layout.fieldNodeCoordinates( Bz, Point{0., 0., 0.}, ix, iy, iz); - Bz(ix, iy, iz) = std::cos(2 * M_PI / 5. * point[0]) - * std::tanh(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + Bz(ix, iy, iz) = std::cos(2 * _PI_ / 5.f * point[0]) + * std::tanh(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -417,7 +451,12 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) for (std::uint32_t iz = psi_p_Z; iz <= pei_p_Z; ++iz) { std::uint32_t index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Jx(ix, iy, iz), ::testing::DoubleNear((expectedJx[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + EXPECT_THAT(Jx(ix, iy, iz), ::testing::DoubleNear((expectedJx[index_]), 1e-12)); + else + { + // todo + } } } } @@ -431,7 +470,12 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) for (std::uint32_t iz = psi_p_Z; iz <= pei_p_Z; ++iz) { std::uint32_t index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Jy(ix, iy, iz), ::testing::DoubleNear((expectedJy[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + EXPECT_THAT(Jy(ix, iy, iz), ::testing::DoubleNear((expectedJy[index_]), 1e-12)); + else + { + // todo + } } } } @@ -445,7 +489,12 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) for (std::uint32_t iz = psi_d_Z; iz <= pei_d_Z; ++iz) { std::uint32_t index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Jz(ix, iy, iz), ::testing::DoubleNear((expectedJz[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + EXPECT_THAT(Jz(ix, iy, iz), ::testing::DoubleNear((expectedJz[index_]), 1e-12)); + else + { + // todo + } } } } diff --git a/tests/core/numerics/faraday/test_main.cpp b/tests/core/numerics/faraday/test_main.cpp index 09637a3d3..8259fd969 100644 --- a/tests/core/numerics/faraday/test_main.cpp +++ b/tests/core/numerics/faraday/test_main.cpp @@ -1,8 +1,6 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" -#include -#include #include "core/data/grid/grid.hpp" @@ -20,6 +18,14 @@ #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" #include "tests/core/data/gridlayout/gridlayout_test.hpp" +#include +#include +#include + + +PHARE::core::floater_t<4> constexpr static _PI_ = std::numbers::pi_v>; + + using namespace PHARE::core; @@ -236,8 +242,8 @@ TEST_F(Faraday1DTest, Faraday1DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ey, Point{0.}, ix); - Ey(ix) = std::cos(2 * M_PI / 5. * point[0]); - Ez(ix) = std::sin(2 * M_PI / 5. * point[0]); + Ey(ix) = std::cos(2 * _PI_ / 5.f * point[0]); + Ez(ix) = std::sin(2 * _PI_ / 5.f * point[0]); } auto gsi_d_X = this->layout.ghostStartIndex(QtyCentering::dual, Direction::X); @@ -247,8 +253,8 @@ TEST_F(Faraday1DTest, Faraday1DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(By, Point{0.}, ix); - By(ix) = std::tanh(point[0] - 5. / 2.); - Bz(ix) = std::tanh(point[0] - 5. / 2.); + By(ix) = std::tanh(point[0] - 5.f / 2.f); + Bz(ix) = std::tanh(point[0] - 5.f / 2.f); } faraday.setLayout(&layout); @@ -259,8 +265,16 @@ TEST_F(Faraday1DTest, Faraday1DCalculatedOk) for (auto ix = psi_d_X; ix <= pei_d_X; ++ix) { - EXPECT_THAT(Bynew(ix), ::testing::DoubleNear((expected_dbydt[ix]), 1e-12)); - EXPECT_THAT(Bznew(ix), ::testing::DoubleNear((expected_dbzdt[ix]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bynew(ix), ::testing::DoubleNear((expected_dbydt[ix]), 1e-12)); + EXPECT_THAT(Bznew(ix), ::testing::DoubleNear((expected_dbzdt[ix]), 1e-12)); + } + else + { + EXPECT_THAT(Bynew(ix), ::testing::FloatNear((expected_dbydt[ix]), 1e-5)); + EXPECT_THAT(Bznew(ix), ::testing::FloatNear((expected_dbzdt[ix]), 1e-5)); + } } } @@ -296,7 +310,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ex, Point{0., 0.}, ix, iy); - Ex(ix, iy) = std::cos(2 * M_PI / 5. * point[0]) * std::sin(2 * M_PI / 6. * point[1]); + Ex(ix, iy) = std::cos(2 * _PI_ / 5.f * point[0]) * std::sin(2 * _PI_ / 6.f * point[1]); } } @@ -306,7 +320,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ey, Point{0., 0.}, ix, iy); - Ey(ix, iy) = std::cos(2 * M_PI / 5. * point[0]) * std::tanh(2 * M_PI / 6. * point[1]); + Ey(ix, iy) = std::cos(2 * _PI_ / 5.f * point[0]) * std::tanh(2 * _PI_ / 6.f * point[1]); } } @@ -316,7 +330,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ez, Point{0., 0.}, ix, iy); - Ez(ix, iy) = std::sin(2 * M_PI / 5. * point[0]) * std::tanh(2 * M_PI / 6. * point[1]); + Ez(ix, iy) = std::sin(2 * _PI_ / 5.f * point[0]) * std::tanh(2 * _PI_ / 6.f * point[1]); } } @@ -326,7 +340,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Bx, Point{0., 0.}, ix, iy); - Bx(ix, iy) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.); + Bx(ix, iy) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f); } } @@ -336,7 +350,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(By, Point{0., 0.}, ix, iy); - By(ix, iy) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.); + By(ix, iy) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f); } } @@ -346,7 +360,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Bz, Point{0., 0.}, ix, iy); - Bz(ix, iy) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.); + Bz(ix, iy) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f); } } @@ -365,7 +379,10 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) for (auto iy = psi_d_Y; iy <= pei_d_Y; ++iy) { auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Bxnew(ix, iy), ::testing::DoubleNear((expected_dbxdt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bxnew(ix, iy), ::testing::DoubleNear((expected_dbxdt[index_]), 1e-12)); + } } } @@ -381,7 +398,14 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) for (auto iy = psi_p_Y; iy <= pei_p_Y; ++iy) { auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Bynew(ix, iy), ::testing::DoubleNear((expected_dbydt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bynew(ix, iy), ::testing::DoubleNear((expected_dbydt[index_]), 1e-12)); + } + else + { + // todo + } } } @@ -392,7 +416,14 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) for (auto iy = psi_d_Y; iy <= pei_d_Y; ++iy) { auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Bznew(ix, iy), ::testing::DoubleNear((expected_dbzdt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bznew(ix, iy), ::testing::DoubleNear((expected_dbzdt[index_]), 1e-12)); + } + else + { + // todo + } } } } @@ -435,9 +466,9 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ex, Point{0., 0., 0.}, ix, iy, iz); - Ex(ix, iy, iz) = std::sin(2 * M_PI / 5. * point[0]) - * std::cos(2 * M_PI / 6. * point[1]) - * std::tanh(2 * M_PI / 12. * point[2]); + Ex(ix, iy, iz) = std::sin(2 * _PI_ / 5.f * point[0]) + * std::cos(2 * _PI_ / 6.f * point[1]) + * std::tanh(2 * _PI_ / 12.f * point[2]); } } } @@ -450,9 +481,9 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ey, Point{0., 0., 0.}, ix, iy, iz); - Ey(ix, iy, iz) = std::tanh(2 * M_PI / 5. * point[0]) - * std::sin(2 * M_PI / 6. * point[1]) - * std::cos(2 * M_PI / 12. * point[2]); + Ey(ix, iy, iz) = std::tanh(2 * _PI_ / 5.f * point[0]) + * std::sin(2 * _PI_ / 6.f * point[1]) + * std::cos(2 * _PI_ / 12.f * point[2]); } } } @@ -465,9 +496,9 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Ez, Point{0., 0., 0.}, ix, iy, iz); - Ez(ix, iy, iz) = std::cos(2 * M_PI / 5. * point[0]) - * std::tanh(2 * M_PI / 6. * point[1]) - * std::sin(2 * M_PI / 12. * point[2]); + Ez(ix, iy, iz) = std::cos(2 * _PI_ / 5.f * point[0]) + * std::tanh(2 * _PI_ / 6.f * point[1]) + * std::sin(2 * _PI_ / 12.f * point[2]); } } } @@ -480,8 +511,8 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Bx, Point{0., 0., 0.}, ix, iy, iz); - Bx(ix, iy, iz) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.) - * std::tanh(point[2] - 12. / 2.); + Bx(ix, iy, iz) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f) + * std::tanh(point[2] - 12.f / 2.f); } } } @@ -494,8 +525,8 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(By, Point{0., 0., 0.}, ix, iy, iz); - By(ix, iy, iz) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.) - * std::tanh(point[2] - 12. / 2.); + By(ix, iy, iz) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f) + * std::tanh(point[2] - 12.f / 2.f); } } } @@ -508,8 +539,8 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) { auto point = this->layout.fieldNodeCoordinates(Bz, Point{0., 0., 0.}, ix, iy, iz); - Bz(ix, iy, iz) = std::tanh(point[0] - 5. / 2.) * std::tanh(point[1] - 6. / 2.) - * std::tanh(point[2] - 12. / 2.); + Bz(ix, iy, iz) = std::tanh(point[0] - 5.f / 2.f) * std::tanh(point[1] - 6.f / 2.f) + * std::tanh(point[2] - 12.f / 2.f); } } } @@ -541,8 +572,11 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) for (auto iz = psi_d_Z; iz <= pei_d_Z; ++iz) { auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Bxnew(ix, iy, iz), - ::testing::DoubleNear((expected_dbxdt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bxnew(ix, iy, iz), + ::testing::DoubleNear((expected_dbxdt[index_]), 1e-12)); + } } } } @@ -556,8 +590,11 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) for (auto iz = psi_d_Z; iz <= pei_d_Z; ++iz) { auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Bynew(ix, iy, iz), - ::testing::DoubleNear((expected_dbydt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bynew(ix, iy, iz), + ::testing::DoubleNear((expected_dbydt[index_]), 1e-12)); + } } } } @@ -571,8 +608,11 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) for (auto iz = psi_p_Z; iz <= pei_p_Z; ++iz) { auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Bznew(ix, iy, iz), - ::testing::DoubleNear((expected_dbzdt[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Bznew(ix, iy, iz), + ::testing::DoubleNear((expected_dbzdt[index_]), 1e-12)); + } } } } diff --git a/tests/core/numerics/interpolator/test_main.cpp b/tests/core/numerics/interpolator/test_main.cpp index 28e2fff14..f94731249 100644 --- a/tests/core/numerics/interpolator/test_main.cpp +++ b/tests/core/numerics/interpolator/test_main.cpp @@ -81,9 +81,10 @@ class AWeighter : public ::testing::Test protected: Weighter weighter; static const int nbr_tests = 100000; - std::array normalizedPositions; - std::array weightsSums; - std::array, nbr_tests> weights_; + std::array, nbr_tests> normalizedPositions; + std::array, nbr_tests> weightsSums; + std::array, nbrPointsSupport(Weighter::interp_order)>, nbr_tests> + weights_; }; @@ -97,8 +98,18 @@ TYPED_TEST_SUITE(AWeighter, Weighters); TYPED_TEST(AWeighter, ComputesWeightThatSumIsOne) { - auto equalsOne = [](double sum) { return std::abs(sum - 1.) < 1e-10; }; - EXPECT_TRUE(std::all_of(std::begin(this->weightsSums), std::end(this->weightsSums), equalsOne)); + if constexpr (std::is_same_v, double>) + { + auto equalsOne = [](auto sum) { return std::abs(sum - 1.) < 1e-10; }; + EXPECT_TRUE( + std::all_of(std::begin(this->weightsSums), std::end(this->weightsSums), equalsOne)); + } + else + { + auto equalsOne = [](auto sum) { return std::abs(sum - 1.f) < 1e-6f; }; + EXPECT_TRUE( + std::all_of(std::begin(this->weightsSums), std::end(this->weightsSums), equalsOne)); + } } @@ -675,7 +686,8 @@ struct ACollectionOfParticles_2d : public ::testing::Test using Grid_t = typename PHARE_TYPES::Grid_t; using UsableVecFieldND = UsableVecField; - GridLayout_t layout{ConstArray(.1), {nx, ny}, ConstArray(0)}; + GridLayout_t layout{ + ConstArray, dim>(.1), {nx, ny}, ConstArray, dim>(0)}; ParticleArray_t particles; Grid_t rho; UsableVecFieldND v; @@ -691,7 +703,7 @@ struct ACollectionOfParticles_2d : public ::testing::Test { auto& part = particles.emplace_back(); part.iCell = {i, j}; - part.delta = ConstArray(.5); + part.delta = ConstArray, dim>(.5); part.weight = 1.; part.v[0] = +2.; part.v[1] = -1.; diff --git a/tests/core/numerics/ion_updater/test_updater.cpp b/tests/core/numerics/ion_updater/test_updater.cpp index 2d939cf7d..9c648cb65 100644 --- a/tests/core/numerics/ion_updater/test_updater.cpp +++ b/tests/core/numerics/ion_updater/test_updater.cpp @@ -11,57 +11,57 @@ using namespace PHARE::core; -using Param = std::vector const&; -using Return = std::shared_ptr>; +using Param = std::vector> const&; +using Return = std::shared_ptr>>; Return density(Param x) { - return std::make_shared>(x.size(), 1); + return std::make_shared>>(x.size(), 1); } Return vx(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } Return vy(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } Return vz(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } Return vthx(Param x) { - return std::make_shared>(x.size(), .1); + return std::make_shared>>(x.size(), .1); } Return vthy(Param x) { - return std::make_shared>(x.size(), .1); + return std::make_shared>>(x.size(), .1); } Return vthz(Param x) { - return std::make_shared>(x.size(), .1); + return std::make_shared>>(x.size(), .1); } Return bx(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } Return by(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } Return bz(Param x) { - return std::make_shared>(x.size(), 0); + return std::make_shared>>(x.size(), 0); } @@ -515,7 +515,7 @@ struct IonUpdaterTest : public ::testing::Test } } // end 1D - } // end pop loop + } // end pop loop PHARE::core::depositParticles(ions, layout, Interpolator{}, PHARE::core::DomainDeposit{}); @@ -574,7 +574,7 @@ struct IonUpdaterTest : public ::testing::Test auto ix1 = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); auto nonZero = [&](auto const& field) { - auto sum = 0.; + floater_t<4> sum = 0; for (auto ix = ix0; ix <= ix1; ++ix) { sum += std::abs(field(ix)); @@ -589,8 +589,8 @@ struct IonUpdaterTest : public ::testing::Test { auto evolution = std::abs(newField(ix) - originalField(ix)); // should check that moments are still compatible with user inputs also - EXPECT_TRUE(evolution > 0.0); - if (evolution <= 0.0) + EXPECT_TRUE(evolution > 0.0f); + if (evolution <= 0.0f) std::cout << "after update : " << newField(ix) << " before update : " << originalField(ix) << " evolution : " << evolution << " ix : " << ix << "\n"; @@ -624,7 +624,7 @@ struct IonUpdaterTest : public ::testing::Test auto check = [&](auto const& density, auto const& function) { std::vector ixes; - std::vector x; + std::vector> x; for (auto ix = ix0; ix < ix1; ++ix) { @@ -642,9 +642,9 @@ struct IonUpdaterTest : public ::testing::Test auto ix = ixes[i]; auto diff = std::abs(density(ix) - functionX[i]); - EXPECT_GE(0.07, diff); + EXPECT_GE(0.07f, diff); - if (diff >= 0.07) + if (diff >= 0.07f) std::cout << "actual : " << density(ix) << " prescribed : " << functionX[i] << " diff : " << diff << " ix : " << ix << "\n"; } diff --git a/tests/core/numerics/ohm/test_main.cpp b/tests/core/numerics/ohm/test_main.cpp index e5d338891..6c35a1fa4 100644 --- a/tests/core/numerics/ohm/test_main.cpp +++ b/tests/core/numerics/ohm/test_main.cpp @@ -68,7 +68,7 @@ struct OhmTest : public ::testing::Test using GridYee = GridLayout>; using UsableVecFieldND = UsableVecField; - using Grid_t = Grid, HybridQuantity::Scalar>; + using Grid_t = Grid>, HybridQuantity::Scalar>; GridYee layout = NDlayout::create(); Grid_t n; @@ -102,23 +102,23 @@ struct OhmTest : public ::testing::Test { auto point = this->layout.fieldNodeCoordinates(n, Point{0.}, ix); - n(ix) = std::cosh(0.5 * point[0]); - Vx(ix) = std::sinh(0.2 * point[0]); - Vy(ix) = std::sinh(0.3 * point[0]); - Vz(ix) = std::sinh(0.4 * point[0]); - P(ix) = std::cosh(0.5 * point[0]); - Bx(ix) = std::cosh(0.2 * point[0]); - Jy(ix) = std::tanh(0.3 * point[0]); - Jz(ix) = std::tanh(0.4 * point[0]); + n(ix) = std::cosh(0.5f * point[0]); + Vx(ix) = std::sinh(0.2f * point[0]); + Vy(ix) = std::sinh(0.3f * point[0]); + Vz(ix) = std::sinh(0.4f * point[0]); + P(ix) = std::cosh(0.5f * point[0]); + Bx(ix) = std::cosh(0.2f * point[0]); + Jy(ix) = std::tanh(0.3f * point[0]); + Jz(ix) = std::tanh(0.4f * point[0]); } for (auto ix = gsi_d_X; ix <= gei_d_X; ++ix) { auto point = this->layout.fieldNodeCoordinates(Bz, Point{0.}, ix); - By(ix) = std::cosh(0.3 * point[0]); - Bz(ix) = std::cosh(0.4 * point[0]); - Jx(ix) = std::tanh(0.2 * point[0]); + By(ix) = std::cosh(0.3f * point[0]); + Bz(ix) = std::cosh(0.4f * point[0]); + Jx(ix) = std::tanh(0.2f * point[0]); } } @@ -140,20 +140,20 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates(n, Point{0., 0.}, ix, iy); - n(ix, iy) = std::cosh(0.5 * point[0]) * std::cosh(0.5 * point[1]); - Vx(ix, iy) = std::sinh(0.2 * point[0]) * std::sinh(0.2 * point[1]); - Vy(ix, iy) = std::sinh(0.3 * point[0]) * std::sinh(0.3 * point[1]); - Vz(ix, iy) = std::sinh(0.4 * point[0]) * std::sinh(0.4 * point[1]); - P(ix, iy) = std::cosh(0.5 * point[0]) * std::cosh(0.5 * point[1]); - Jz(ix, iy) = std::tanh(0.4 * point[0]) * std::tanh(0.4 * point[1]); + n(ix, iy) = std::cosh(0.5f * point[0]) * std::cosh(0.5f * point[1]); + Vx(ix, iy) = std::sinh(0.2f * point[0]) * std::sinh(0.2f * point[1]); + Vy(ix, iy) = std::sinh(0.3f * point[0]) * std::sinh(0.3f * point[1]); + Vz(ix, iy) = std::sinh(0.4f * point[0]) * std::sinh(0.4f * point[1]); + P(ix, iy) = std::cosh(0.5f * point[0]) * std::cosh(0.5f * point[1]); + Jz(ix, iy) = std::tanh(0.4f * point[0]) * std::tanh(0.4f * point[1]); } for (auto iy = gsi_d_Y; iy <= gei_d_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(Bx, Point{0., 0.}, ix, iy); - Bx(ix, iy) = std::cosh(0.2 * point[0]) * std::cosh(0.2 * point[1]); - Jy(ix, iy) = std::tanh(0.3 * point[0]) * std::tanh(0.3 * point[1]); + Bx(ix, iy) = std::cosh(0.2f * point[0]) * std::cosh(0.2f * point[1]); + Jy(ix, iy) = std::tanh(0.3f * point[0]) * std::tanh(0.3f * point[1]); } } for (auto ix = gsi_d_X; ix <= gei_d_X; ++ix) @@ -163,15 +163,15 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates(Jx, Point{0., 0.}, ix, iy); - By(ix, iy) = std::cosh(0.3 * point[0]) * std::cosh(0.3 * point[1]); - Jx(ix, iy) = std::tanh(0.2 * point[0]) * std::tanh(0.2 * point[1]); + By(ix, iy) = std::cosh(0.3f * point[0]) * std::cosh(0.3f * point[1]); + Jx(ix, iy) = std::tanh(0.2f * point[0]) * std::tanh(0.2f * point[1]); } for (auto iy = gsi_d_Y; iy <= gei_d_Y; ++iy) { auto point = this->layout.fieldNodeCoordinates(Bz, Point{0., 0.}, ix, iy); - Bz(ix, iy) = std::cosh(0.4 * point[0]) * std::cosh(0.4 * point[1]); + Bz(ix, iy) = std::cosh(0.4f * point[0]) * std::cosh(0.4f * point[1]); } } } @@ -200,24 +200,24 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates( n, Point{0., 0., 0.}, ix, iy, iz); - n(ix, iy, iz) = std::cosh(0.5 * point[0]) * std::cosh(0.5 * point[1]) - * std::cosh(0.5 * point[2]); - Vx(ix, iy, iz) = std::sinh(0.2 * point[0]) * std::sinh(0.2 * point[1]) - * std::sinh(0.2 * point[2]); - Vy(ix, iy, iz) = std::sinh(0.3 * point[0]) * std::sinh(0.3 * point[1]) - * std::sinh(0.3 * point[2]); - Vz(ix, iy, iz) = std::sinh(0.4 * point[0]) * std::sinh(0.4 * point[1]) - * std::sinh(0.4 * point[2]); - P(ix, iy, iz) = std::cosh(0.5 * point[0]) * std::cosh(0.5 * point[1]) - * std::cosh(0.5 * point[2]); + n(ix, iy, iz) = std::cosh(0.5f * point[0]) * std::cosh(0.5f * point[1]) + * std::cosh(0.5f * point[2]); + Vx(ix, iy, iz) = std::sinh(0.2f * point[0]) * std::sinh(0.2f * point[1]) + * std::sinh(0.2f * point[2]); + Vy(ix, iy, iz) = std::sinh(0.3f * point[0]) * std::sinh(0.3f * point[1]) + * std::sinh(0.3f * point[2]); + Vz(ix, iy, iz) = std::sinh(0.4f * point[0]) * std::sinh(0.4f * point[1]) + * std::sinh(0.4f * point[2]); + P(ix, iy, iz) = std::cosh(0.5f * point[0]) * std::cosh(0.5f * point[1]) + * std::cosh(0.5f * point[2]); } for (auto iz = gsi_d_Z; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates( Jz, Point{0., 0., 0.}, ix, iy, iz); - Jz(ix, iy, iz) = std::tanh(0.4 * point[0]) * std::tanh(0.4 * point[1]) - * std::tanh(0.4 * point[2]); + Jz(ix, iy, iz) = std::tanh(0.4f * point[0]) * std::tanh(0.4f * point[1]) + * std::tanh(0.4f * point[2]); } } for (auto iy = gsi_d_Y; iy <= gei_d_Y; ++iy) @@ -227,16 +227,16 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates( Jy, Point{0., 0., 0.}, ix, iy, iz); - Jy(ix, iy, iz) = std::tanh(0.3 * point[0]) * std::tanh(0.3 * point[1]) - * std::tanh(0.3 * point[2]); + Jy(ix, iy, iz) = std::tanh(0.3f * point[0]) * std::tanh(0.3f * point[1]) + * std::tanh(0.3f * point[2]); } for (auto iz = gsi_d_Z; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates( Bx, Point{0., 0., 0.}, ix, iy, iz); - Bx(ix, iy, iz) = std::cosh(0.2 * point[0]) * std::cosh(0.2 * point[1]) - * std::cosh(0.2 * point[2]); + Bx(ix, iy, iz) = std::cosh(0.2f * point[0]) * std::cosh(0.2f * point[1]) + * std::cosh(0.2f * point[2]); } } } @@ -249,16 +249,16 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates( Jx, Point{0., 0., 0.}, ix, iy, iz); - Jx(ix, iy, iz) = std::tanh(0.2 * point[0]) * std::tanh(0.2 * point[1]) - * std::tanh(0.2 * point[2]); + Jx(ix, iy, iz) = std::tanh(0.2f * point[0]) * std::tanh(0.2f * point[1]) + * std::tanh(0.2f * point[2]); } for (auto iz = gsi_d_Z; iz <= gei_d_Z; ++iz) { auto point = this->layout.fieldNodeCoordinates( By, Point{0., 0., 0.}, ix, iy, iz); - By(ix, iy, iz) = std::cosh(0.3 * point[0]) * std::cosh(0.3 * point[1]) - * std::cosh(0.3 * point[2]); + By(ix, iy, iz) = std::cosh(0.3f * point[0]) * std::cosh(0.3f * point[1]) + * std::cosh(0.3f * point[2]); } } for (auto iy = gsi_d_Y; iy <= gei_d_Y; ++iy) @@ -268,8 +268,8 @@ struct OhmTest : public ::testing::Test auto point = this->layout.fieldNodeCoordinates( Bz, Point{0., 0., 0.}, ix, iy, iz); - Bz(ix, iy, iz) = std::cosh(0.4 * point[0]) * std::cosh(0.4 * point[1]) - * std::cosh(0.4 * point[2]); + Bz(ix, iy, iz) = std::cosh(0.4f * point[0]) * std::cosh(0.4f * point[1]) + * std::cosh(0.4f * point[2]); } } } @@ -321,7 +321,7 @@ TYPED_TEST(OhmTest, ShouldBeGivenAGridLayoutPointerToBeOperational) } -std::vector read(std::string filename) +auto read(std::string filename) { std::ifstream readFile(filename); assert(readFile.is_open()); @@ -329,7 +329,16 @@ std::vector read(std::string filename) std::copy(std::istream_iterator(readFile), std::istream_iterator(), std::back_inserter(x)); - return x; + + if constexpr (std::is_same_v, float>) + { + std::vector> fx(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) + fx[i] = x[i]; + return fx; + } + else + return x; } @@ -368,7 +377,10 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) for (auto ix = psi_X; ix <= pei_X; ++ix) { - EXPECT_THAT(Exnew(ix), ::testing::DoubleNear((expected_ohmX[ix]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Exnew(ix), ::testing::DoubleNear((expected_ohmX[ix]), 1e-12)); + } } psi_X = this->layout.physicalStartIndex(Eynew, Direction::X); @@ -376,7 +388,10 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) for (auto ix = psi_X; ix <= pei_X; ++ix) { - EXPECT_THAT(Eynew(ix), ::testing::DoubleNear((expected_ohmY[ix]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eynew(ix), ::testing::DoubleNear((expected_ohmY[ix]), 1e-12)); + } } psi_X = this->layout.physicalStartIndex(Eznew, Direction::X); @@ -384,7 +399,10 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) for (auto ix = psi_X; ix <= pei_X; ++ix) { - EXPECT_THAT(Eznew(ix), ::testing::DoubleNear((expected_ohmZ[ix]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eznew(ix), ::testing::DoubleNear((expected_ohmZ[ix]), 1e-12)); + } } } @@ -401,7 +419,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ex); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Exnew(ix, iy), ::testing::DoubleNear((expected_ohmX[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Exnew(ix, iy), + ::testing::DoubleNear((expected_ohmX[index_]), 1e-12)); + } } } @@ -416,7 +438,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ey); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Eynew(ix, iy), ::testing::DoubleNear((expected_ohmY[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eynew(ix, iy), + ::testing::DoubleNear((expected_ohmY[index_]), 1e-12)); + } } } @@ -431,7 +457,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ez); auto index_ = ix * nPts_[1] + iy; - EXPECT_THAT(Eznew(ix, iy), ::testing::DoubleNear((expected_ohmZ[index_]), 1e-12)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eznew(ix, iy), + ::testing::DoubleNear((expected_ohmZ[index_]), 1e-12)); + } } } } @@ -453,8 +483,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ex); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Exnew(ix, iy, iz), - ::testing::DoubleNear((expected_ohmX[index_]), 1e-10)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Exnew(ix, iy, iz), + ::testing::DoubleNear((expected_ohmX[index_]), 1e-10)); + } } } } @@ -474,8 +507,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ey); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Eynew(ix, iy, iz), - ::testing::DoubleNear((expected_ohmY[index_]), 1e-10)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eynew(ix, iy, iz), + ::testing::DoubleNear((expected_ohmY[index_]), 1e-10)); + } } } } @@ -495,8 +531,11 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { auto nPts_ = this->layout.allocSize(HybridQuantity::Scalar::Ez); auto index_ = ix * nPts_[1] * nPts_[2] + iy * nPts_[2] + iz; - EXPECT_THAT(Eznew(ix, iy, iz), - ::testing::DoubleNear((expected_ohmZ[index_]), 1e-10)); + if constexpr (std::is_same_v, double>) + { + EXPECT_THAT(Eznew(ix, iy, iz), + ::testing::DoubleNear((expected_ohmZ[index_]), 1e-10)); + } } } } diff --git a/tests/core/numerics/pusher/test_pusher.cpp b/tests/core/numerics/pusher/test_pusher.cpp index 8bad3a9d1..f0290bc02 100644 --- a/tests/core/numerics/pusher/test_pusher.cpp +++ b/tests/core/numerics/pusher/test_pusher.cpp @@ -2,6 +2,7 @@ #include "gtest/gtest.h" #include +#include #include #include #include @@ -125,6 +126,8 @@ struct DummyLayout template class APusher : public ::testing::Test { + using ParticleArray_t = ParticleArray; + public: using Pusher_ = BorisPusher>, Electromag, Interpolator, BoundaryCondition, DummyLayout>; @@ -141,10 +144,9 @@ class APusher : public ::testing::Test , nt{static_cast((tend - tstart) / dt + 1)} { - particlesIn.emplace_back( - Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); - particlesOut.emplace_back( - Particle{1., 1., ConstArray(5), ConstArray(0.), {0., 10., 0.}}); + particlesIn.emplace_back(Particle{ + 1., 1., ConstArray(5), ConstArray, dim>(0.), {0., 10., 0.}}); + particlesOut.emplace_back(particlesIn.back()); dxyz.fill(0.05); for (std::size_t i = 0; i < dim; i++) actual[i].resize(nt, 0.05); @@ -168,7 +170,7 @@ class APusher : public ::testing::Test // BoundaryCondition bc; std::array, dim> actual; - std::array dxyz; + std::array, dim> dxyz; }; @@ -184,9 +186,12 @@ TEST_F(APusher3D, trajectoryIsOk) for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; - actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1]; - actual[2][i] = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * dxyz[2]; + actual[0][i] + = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * floater_t<0>(dxyz[0]); + actual[1][i] + = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * floater_t<0>(dxyz[1]); + actual[2][i] + = (particlesOut[0].iCell[2] + particlesOut[0].delta[2]) * floater_t<0>(dxyz[2]); pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); @@ -209,8 +214,10 @@ TEST_F(APusher2D, trajectoryIsOk) for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; - actual[1][i] = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * dxyz[1]; + actual[0][i] + = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * floater_t<0>(dxyz[0]); + actual[1][i] + = (particlesOut[0].iCell[1] + particlesOut[0].delta[1]) * floater_t<0>(dxyz[1]); pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); @@ -231,7 +238,8 @@ TEST_F(APusher1D, trajectoryIsOk) for (decltype(nt) i = 0; i < nt; ++i) { - actual[0][i] = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * dxyz[0]; + actual[0][i] + = (particlesOut[0].iCell[0] + particlesOut[0].delta[0]) * floater_t<0>(dxyz[0]); pusher->move(rangeIn, rangeOut, em, mass, interpolator, layout, selector, selector); @@ -249,11 +257,13 @@ TEST_F(APusher1D, trajectoryIsOk) // and those that stay. class APusherWithLeavingParticles : public ::testing::Test { + using ParticleArray_t = ParticleArray<1>; + public: APusherWithLeavingParticles() : pusher{std::make_unique< - BorisPusher<1, IndexRange>, Electromag, Interpolator, - BoundaryCondition<1, 1>, DummyLayout<1>>>()} + BorisPusher<1, IndexRange>, Electromag, Interpolator, + BoundaryCondition<1, 1>, DummyLayout<1>>>()} , mass{1} , dt{0.001} , tstart{0} @@ -268,7 +278,7 @@ class APusherWithLeavingParticles : public ::testing::Test std::random_device rd; std::mt19937 gen(rd()); std::uniform_int_distribution<> dis(0, 9); - std::uniform_real_distribution delta(0, 1); + std::uniform_real_distribution> delta(0, 1); for (std::size_t iPart = 0; iPart < 1000; ++iPart) { @@ -290,13 +300,13 @@ class APusherWithLeavingParticles : public ::testing::Test std::size_t nt; Electromag em; Interpolator interpolator; - double dx = 0.1; + floater_t<0> dx = 0.1; Box domain; Box cells; BoundaryCondition<1, 1> bc; - ParticleArray<1> particlesIn; - ParticleArray<1> particlesOut1; - ParticleArray<1> particlesOut2; + ParticleArray_t particlesIn; + ParticleArray_t particlesOut1; + ParticleArray_t particlesOut2; }; diff --git a/tests/initializer/init_functions.hpp b/tests/initializer/init_functions.hpp index 9a7ad3c6b..317975bb2 100644 --- a/tests/initializer/init_functions.hpp +++ b/tests/initializer/init_functions.hpp @@ -6,59 +6,63 @@ #include "core/utilities/span.hpp" + namespace PHARE::initializer::test_fn::func_1d { -using Param = std::vector const&; -using Return = std::shared_ptr>; + +using namespace PHARE::core; + +using Param = std::vector> const&; +using Return = std::shared_ptr>>; Return density(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vx(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vy(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vz(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthx(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthy(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthz(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return bx(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return by(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return bz(Param x) { - return std::make_shared>(x); + return std::make_shared>>(x); } @@ -67,57 +71,60 @@ Return bz(Param x) namespace PHARE::initializer::test_fn::func_2d { -using Param = std::vector const&; -using Return = std::shared_ptr>; + +using namespace PHARE::core; + +using Param = std::vector> const&; +using Return = std::shared_ptr>>; Return density(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vx(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vy(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vz(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthx(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthy(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return vthz(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return bx(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return by(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } Return bz(Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); } @@ -127,22 +134,24 @@ Return bz(Param x, Param /*y*/) template auto makeSharedPtr() { - using Param = std::vector const&; + using namespace PHARE::core; + + using Param = std::vector> const&; if constexpr (dim == 1) { - return [](Param x) { return std::make_shared>(x); }; + return [](Param x) { return std::make_shared>>(x); }; } else if constexpr (dim == 2) { return [](Param x, Param /*y*/) { - return std::make_shared>(x); + return std::make_shared>>(x); }; } else if constexpr (dim == 3) { return [](Param x, Param /*y*/, Param /*z*/) { - return std::make_shared>(x); + return std::make_shared>>(x); }; } } diff --git a/tests/initializer/test_initializer.cpp b/tests/initializer/test_initializer.cpp index 2215eab68..5c4e4b29b 100644 --- a/tests/initializer/test_initializer.cpp +++ b/tests/initializer/test_initializer.cpp @@ -17,7 +17,7 @@ #include "core/data/particles/particle_array.hpp" - +using namespace PHARE::core; using namespace PHARE::initializer; using GridLayoutT = PHARE::core::GridLayout>; @@ -57,12 +57,12 @@ TEST(APythonDataProvider, providesAValidTree) auto simulationName = input["simulation"]["name"].to(); auto dim = input["simulation"]["dimension"].to(); auto interp_order = input["simulation"]["interp_order"].to(); - auto dt = input["simulation"]["time_step"].to(); + auto dt = input["simulation"]["time_step"].to>(); auto layout = input["simulation"]["grid"]["layout_type"].to(); auto nx = input["simulation"]["grid"]["nbr_cells"]["x"].to(); - auto dx = input["simulation"]["grid"]["meshsize"]["x"].to(); - auto origin = input["simulation"]["grid"]["origin"]["x"].to(); + auto dx = input["simulation"]["grid"]["meshsize"]["x"].to>(); + auto origin = input["simulation"]["grid"]["origin"]["x"].to>(); auto pusherName = input["simulation"]["algo"]["ion_updater"]["pusher"]["name"].to(); @@ -71,7 +71,7 @@ TEST(APythonDataProvider, providesAValidTree) auto nbrPopulations = input["simulation"]["ions"]["nbrPopulations"].to(); auto& pop0 = input["simulation"]["ions"]["pop0"]; auto pop0Name = pop0["name"].to(); - auto pop0Mass = pop0["mass"].to(); + auto pop0Mass = pop0["mass"].to>(); auto& pop0ParticleInitializer = pop0["particle_initializer"]; auto pop0ParticleInitializerName = pop0ParticleInitializer["name"].to(); auto pop0density = pop0ParticleInitializer["density"].to>(); @@ -82,7 +82,7 @@ TEST(APythonDataProvider, providesAValidTree) auto vth0y = pop0ParticleInitializer["thermal_velocity_y"].to>(); auto vth0z = pop0ParticleInitializer["thermal_velocity_z"].to>(); auto pop0NbrPartPerCell = pop0ParticleInitializer["nbr_part_per_cell"].to(); - auto pop0Charge = pop0ParticleInitializer["charge"].to(); + auto pop0Charge = pop0ParticleInitializer["charge"].to>(); auto pop0Basis = pop0ParticleInitializer["basis"].to(); @@ -91,28 +91,28 @@ TEST(APythonDataProvider, providesAValidTree) EXPECT_EQ(1, dim); EXPECT_EQ(65, nx); - EXPECT_DOUBLE_EQ(1. / 65., dx); - EXPECT_DOUBLE_EQ(0.001, dt); - EXPECT_DOUBLE_EQ(0., origin); + EXPECT_FLOAT_EQ(1. / 65., dx); + EXPECT_FLOAT_EQ(0.001, dt); + EXPECT_FLOAT_EQ(0., origin); EXPECT_EQ("yee", layout); EXPECT_EQ("modified_boris", pusherName); - std::vector input_2 = {2}; + std::vector> input_2 = {2}; EXPECT_EQ(2, nbrPopulations); EXPECT_EQ("protons", pop0Name); - EXPECT_DOUBLE_EQ(1., pop0Mass); + EXPECT_FLOAT_EQ(1., pop0Mass); EXPECT_EQ("maxwellian", pop0ParticleInitializerName); - EXPECT_DOUBLE_EQ(2., (*pop0density(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*bulk0x(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*bulk0y(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*bulk0z(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*vth0x(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*vth0y(input_2))[0]); - EXPECT_DOUBLE_EQ(1., (*vth0z(input_2))[0]); + EXPECT_FLOAT_EQ(2., (*pop0density(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*bulk0x(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*bulk0y(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*bulk0z(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*vth0x(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*vth0y(input_2))[0]); + EXPECT_FLOAT_EQ(1., (*vth0z(input_2))[0]); EXPECT_EQ(100, pop0NbrPartPerCell); - EXPECT_DOUBLE_EQ(1., pop0Charge); + EXPECT_FLOAT_EQ(1., pop0Charge); EXPECT_EQ("cartesian", pop0Basis); PHAREDictHandler::INSTANCE().stop(); diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 569f2c873..c21f228d6 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -289,7 +289,9 @@ def base_test_overlaped_fields_are_equal(self, datahier, coarsest_time): # seems correct considering ghosts are filled with schedules # involving linear/spatial interpolations and so on where # rounding errors may occur.... setting atol to 5.5e-15 - assert_fp_any_all_close(slice1, slice2, atol=5.5e-15, rtol=0) + assert_fp_any_all_close( + slice1, slice2, atol=5.5e-15, atol_fp32=8e-07, rtol=0 + ) checks += 1 except AssertionError as e: print("AssertionError", pd1.name, e) diff --git a/tools/bench/core/bench.hpp b/tools/bench/core/bench.hpp index c46df93a6..472e3d13c 100644 --- a/tools/bench/core/bench.hpp +++ b/tools/bench/core/bench.hpp @@ -21,7 +21,7 @@ PHARE::core::Particle particle(int icell = 15) return {/*.weight = */ 0, /*.charge = */ 1, /*.iCell = */ PHARE::core::ConstArray(icell), - /*.delta = */ PHARE::core::ConstArray(.5), + /*.delta = */ PHARE::core::ConstArray, dim>(.5), /*.v = */ {{.00001, .00001, .00001}}}; }