Skip to content

Commit

Permalink
fix for PHAREHUB#176 - incorrect interpolation in 2/3D
Browse files Browse the repository at this point in the history
  • Loading branch information
@PhilipDeegan
PhilipDeegan committed Jul 16, 2020
1 parent 5cf58ec commit d83a19d
Show file tree
Hide file tree
Showing 3 changed files with 97 additions and 45 deletions.
50 changes: 26 additions & 24 deletions src/core/numerics/interpolator/interpolator.h
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#ifndef PHARE_CORE_NUMERICS_INTERPOLATOR_INTERPOLATOR_H
#define PHARE_CORE_NUMERICS_INTERPOLATOR_INTERPOLATOR_H



#include <array>
#include <cstddef>

Expand Down Expand Up @@ -360,17 +362,17 @@ namespace core



auto const& yDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][0];
auto const& yDenWeights = weights[static_cast<int>(densityCentering[0])][0];
auto const& yDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][1];
auto const& yDenWeights = weights[static_cast<int>(densityCentering[0])][1];

auto const& yXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][1])][0];
auto const& yXFluxWeights = weights[static_cast<int>(fluxCentering[0][1])][0];
auto const& yXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][1])][1];
auto const& yXFluxWeights = weights[static_cast<int>(fluxCentering[0][1])][1];

auto const& yYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][1])][0];
auto const& yYFluxWeights = weights[static_cast<int>(fluxCentering[1][1])][0];
auto const& yYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][1])][1];
auto const& yYFluxWeights = weights[static_cast<int>(fluxCentering[1][1])][1];

auto const& yZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][1])][0];
auto const& yZFluxWeights = weights[static_cast<int>(fluxCentering[2][1])][0];
auto const& yZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][1])][1];
auto const& yZFluxWeights = weights[static_cast<int>(fluxCentering[2][1])][1];



Expand Down Expand Up @@ -441,32 +443,32 @@ namespace core



auto const& yDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][0];
auto const& yDenWeights = weights[static_cast<int>(densityCentering[0])][0];
auto const& yDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][1];
auto const& yDenWeights = weights[static_cast<int>(densityCentering[0])][1];

auto const& yXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][1])][0];
auto const& yXFluxWeights = weights[static_cast<int>(fluxCentering[0][1])][0];
auto const& yXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][1])][1];
auto const& yXFluxWeights = weights[static_cast<int>(fluxCentering[0][1])][1];

auto const& yYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][1])][0];
auto const& yYFluxWeights = weights[static_cast<int>(fluxCentering[1][1])][0];
auto const& yYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][1])][1];
auto const& yYFluxWeights = weights[static_cast<int>(fluxCentering[1][1])][1];

auto const& yZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][1])][0];
auto const& yZFluxWeights = weights[static_cast<int>(fluxCentering[2][1])][0];
auto const& yZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][1])][1];
auto const& yZFluxWeights = weights[static_cast<int>(fluxCentering[2][1])][1];




auto const& zDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][0];
auto const& zDenWeights = weights[static_cast<int>(densityCentering[0])][0];
auto const& zDenStartIndex = startIndex[static_cast<int>(densityCentering[0])][2];
auto const& zDenWeights = weights[static_cast<int>(densityCentering[0])][2];

auto const& zXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][2])][0];
auto const& zXFluxWeights = weights[static_cast<int>(fluxCentering[0][2])][0];
auto const& zXFluxStartIndex = startIndex[static_cast<int>(fluxCentering[0][2])][2];
auto const& zXFluxWeights = weights[static_cast<int>(fluxCentering[0][2])][2];

auto const& zYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][2])][0];
auto const& zYFluxWeights = weights[static_cast<int>(fluxCentering[1][2])][0];
auto const& zYFluxStartIndex = startIndex[static_cast<int>(fluxCentering[1][2])][2];
auto const& zYFluxWeights = weights[static_cast<int>(fluxCentering[1][2])][2];

auto const& zZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][2])][0];
auto const& zZFluxWeights = weights[static_cast<int>(fluxCentering[2][2])][0];
auto const& zZFluxStartIndex = startIndex[static_cast<int>(fluxCentering[2][2])][2];
auto const& zZFluxWeights = weights[static_cast<int>(fluxCentering[2][2])][2];

auto const partRho = particle.weight * coef / cellVolume;
auto const xPartFlux = particle.v[0] * particle.weight * coef / cellVolume;
Expand Down
9 changes: 9 additions & 0 deletions src/core/utilities/types.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,15 @@ namespace core
return sized_array<To_Size>(arr);
}

template<typename Type, size_t size>
constexpr decltype(auto) ConstArray(Type val = 0)
{
std::array<Type, size> arr{};
for (uint8_t i = 0; i < size; i++)
arr[i] = val;
return arr;
}


} // namespace core
} // namespace PHARE
Expand Down
83 changes: 62 additions & 21 deletions tests/core/numerics/interpolator/test_main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@


#include "gmock/gmock.h"
#include "gtest/gtest.h"

Expand Down Expand Up @@ -525,7 +527,7 @@ TYPED_TEST(A3DInterpolator, canComputeAllEMfieldsAtParticle)


template<typename Interpolator>
class ACollectionOfParticles : public ::testing::Test
class ACollectionOfParticles_1d : public ::testing::Test
{
public:
static constexpr uint32_t nx = 40;
Expand All @@ -543,7 +545,7 @@ class ACollectionOfParticles : public ::testing::Test



ACollectionOfParticles()
ACollectionOfParticles_1d()
: part{}
, particles{}
, rho{"field", HybridQuantity::Scalar::rho, nx}
Expand Down Expand Up @@ -669,38 +671,77 @@ class ACollectionOfParticles : public ::testing::Test
protected:
Interpolator interpolator;
};
TYPED_TEST_SUITE_P(ACollectionOfParticles_1d);



TYPED_TEST_SUITE_P(ACollectionOfParticles);


TYPED_TEST_P(ACollectionOfParticles, DepositCorrectlyTheirWeight)
TYPED_TEST_P(ACollectionOfParticles_1d, DepositCorrectlyTheirWeight_1d)
{
EXPECT_DOUBLE_EQ(this->rho(25u), 1.0);
EXPECT_DOUBLE_EQ(this->vx(25u), 2.0);
EXPECT_DOUBLE_EQ(this->vy(25u), -1.0);
EXPECT_DOUBLE_EQ(this->vz(25u), 1.0);
EXPECT_DOUBLE_EQ(this->rho(25), 1.0);
EXPECT_DOUBLE_EQ(this->vx(25), 2.0);
EXPECT_DOUBLE_EQ(this->vy(25), -1.0);
EXPECT_DOUBLE_EQ(this->vz(25), 1.0);
}
REGISTER_TYPED_TEST_SUITE_P(ACollectionOfParticles_1d, DepositCorrectlyTheirWeight_1d);

using MyTypes = ::testing::Types<Interpolator<1, 1>, Interpolator<1, 2>, Interpolator<1, 3>>;
INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_1d, MyTypes);


template<typename Interpolator>
struct ACollectionOfParticles_2d : public ::testing::Test
{
static constexpr size_t dim = 2;
static constexpr uint32_t nx = 15, ny = 15;
static constexpr int start = 0, end = 5;

ACollectionOfParticles_2d()
: rho{"field", HybridQuantity::Scalar::rho, nx, ny}
, vx{"v_x", HybridQuantity::Scalar::Vx, nx, ny}
, vy{"v_y", HybridQuantity::Scalar::Vy, nx, ny}
, vz{"v_z", HybridQuantity::Scalar::Vz, nx, ny}
, v{"v", HybridQuantity::Vector::V}
{
v.setBuffer("v_x", &vx);
v.setBuffer("v_y", &vy);
v.setBuffer("v_z", &vz);

REGISTER_TYPED_TEST_SUITE_P(ACollectionOfParticles, DepositCorrectlyTheirWeight);
for (int i = start; i < end; i++)
for (int j = start; j < end; j++)
{
auto& part = particles.emplace_back();
part.iCell = {i, j};
part.delta = ConstArray<float, dim>(.5);
part.weight = .1 * layout.meshSize()[0];
part.v[0] = +2.;
part.v[1] = -1.;
part.v[2] = +1.;
}

interpolator(std::begin(particles), std::end(particles), rho, v, layout);
}

GridLayout<GridLayoutImplYee<dim, Interpolator::interp_order>> layout{
ConstArray<double, dim>(.1), {nx, ny}, ConstArray<double, dim>(0)};

/*using GridLayoutYee2DO1 = GridLayout<GridLayoutImplYee<2,1>>;
using GridLayoutYee2DO2 = GridLayout<GridLayoutImplYee<2,2>>;
using GridLayoutYee2DO3 = GridLayout<GridLayoutImplYee<2,3>>;
using GridLayoutYee3DO1 = GridLayout<GridLayoutImplYee<3,1>>;
using GridLayoutYee3DO2 = GridLayout<GridLayoutImplYee<3,2>>;
using GridLayoutYee3DO3 = GridLayout<GridLayoutImplYee<3,3>>;*/
ParticleArray<dim> particles;
Field<NdArrayVector<dim>, typename HybridQuantity::Scalar> rho, vx, vy, vz;
VecField<NdArrayVector<dim>, HybridQuantity> v;
Interpolator interpolator;
};
TYPED_TEST_SUITE_P(ACollectionOfParticles_2d);


TYPED_TEST_P(ACollectionOfParticles_2d, DepositCorrectlyTheirWeight_2d)
{
EXPECT_DOUBLE_EQ(this->rho(7, 7), 1.0);
EXPECT_DOUBLE_EQ(this->vx(7, 7), 2.0);
EXPECT_DOUBLE_EQ(this->vy(7, 7), -1.0);
EXPECT_DOUBLE_EQ(this->vz(7, 7), 1.0);
}
REGISTER_TYPED_TEST_SUITE_P(ACollectionOfParticles_2d, DepositCorrectlyTheirWeight_2d);

using MyTypes = ::testing::Types<Interpolator<1, 1>, Interpolator<1, 2>, Interpolator<1, 3>>;
INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles, MyTypes);

using My2dTypes = ::testing::Types<Interpolator<2, 1>, Interpolator<2, 2>, Interpolator<2, 3>>;
INSTANTIATE_TYPED_TEST_SUITE_P(testInterpolator, ACollectionOfParticles_2d, My2dTypes);


int main(int argc, char** argv)
Expand Down

0 comments on commit d83a19d

Please sign in to comment.