Skip to content

Commit

Permalink
variable nu
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolasaunai
nicolasaunai committed Sep 19, 2024
1 parent ce8a371 commit 9a12fd0
Show file tree
Hide file tree
Showing 4 changed files with 175 additions and 8 deletions.
5 changes: 5 additions & 0 deletions src/core/data/grid/gridlayout.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1073,6 +1073,7 @@ namespace core
*/
NO_DISCARD auto static constexpr JzToMoments() { return GridLayoutImpl::JzToMoments(); }

NO_DISCARD auto static constexpr BxToEx() { return GridLayoutImpl::BxToEx(); }

/**
* @brief ByToEx return the indexes and associated coef to compute the linear
Expand All @@ -1088,13 +1089,15 @@ namespace core
NO_DISCARD auto static constexpr BzToEx() { return GridLayoutImpl::BzToEx(); }



/**
* @brief BxToEy return the indexes and associated coef to compute the linear
* interpolation necessary to project Bx onto Ey.
*/
NO_DISCARD auto static constexpr BxToEy() { return GridLayoutImpl::BxToEy(); }


NO_DISCARD auto static constexpr ByToEy() { return GridLayoutImpl::ByToEy(); }

/**
* @brief BzToEy return the indexes and associated coef to compute the linear
Expand All @@ -1118,6 +1121,8 @@ namespace core
*/
NO_DISCARD auto static constexpr ByToEz() { return GridLayoutImpl::ByToEz(); }

NO_DISCARD auto static constexpr BzToEz() { return GridLayoutImpl::BzToEz(); }



/**
Expand Down
122 changes: 122 additions & 0 deletions src/core/data/grid/gridlayoutimplyee.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -685,6 +685,47 @@ namespace core
}
}

NO_DISCARD auto static constexpr BxToEx()
{
// Bx is primal dual dual
// Ex is dual primal primal
// operation is pdd to dpp
[[maybe_unused]] auto constexpr p2dShift = primalToDual();
[[maybe_unused]] auto constexpr d2pShift = dualToPrimal();

if constexpr (dimension == 1)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{p2dShift}, 0.5};
return std::array<WeightPoint<dimension>, 2>{P1, P2};
}
else if constexpr (dimension == 2)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.25};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, d2pShift}, 0.25};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{p2dShift, 0}, 0.25};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{p2dShift, d2pShift},
0.25};
return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
}
else if constexpr (dimension == 3)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.125};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{0, d2pShift, 0}, 0.125};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{p2dShift, 0, 0}, 0.125};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{p2dShift, d2pShift, 0},
0.125};

constexpr WeightPoint<dimension> P5{Point<int, dimension>{0, 0, d2pShift}, 0.125};
constexpr WeightPoint<dimension> P6{Point<int, dimension>{0, d2pShift, d2pShift},
0.125};
constexpr WeightPoint<dimension> P7{Point<int, dimension>{p2dShift, 0, d2pShift},
0.125};
constexpr WeightPoint<dimension> P8{
Point<int, dimension>{p2dShift, d2pShift, d2pShift}, 0.125};
return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
}
}


NO_DISCARD auto static constexpr ByToEx()
Expand Down Expand Up @@ -744,6 +785,47 @@ namespace core
}


NO_DISCARD auto static constexpr BzToEz()
{
// Bz is dual dual primal
// Ez is primal primal dual
// operation is thus ddp to ppd
auto constexpr p2dShift = primalToDual();
auto constexpr d2pShift = dualToPrimal();

if constexpr (dimension == 1)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift}, 0.5};
return std::array<WeightPoint<dimension>, 2>{P1, P2};
}

else if constexpr (dimension == 2)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.25};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift, 0}, 0.25};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, d2pShift}, 0.25};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, d2pShift},
0.25};
return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
}
else if constexpr (dimension == 3)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.25};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift, 0, 0}, 0.25};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, d2pShift, 0}, 0.25};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, d2pShift, 0},
0.25};
constexpr WeightPoint<dimension> P5{Point<int, dimension>{0, 0, p2dShift}, 0.25};
constexpr WeightPoint<dimension> P6{Point<int, dimension>{d2pShift, 0, p2dShift},
0.25};
constexpr WeightPoint<dimension> P7{Point<int, dimension>{0, d2pShift, p2dShift},
0.25};
constexpr WeightPoint<dimension> P8{
Point<int, dimension>{d2pShift, d2pShift, p2dShift}, 0.25};
return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
}
}


NO_DISCARD auto static constexpr ByToEz()
Expand Down Expand Up @@ -838,6 +920,46 @@ namespace core
}


NO_DISCARD auto static constexpr ByToEy()
{
// By is dual primal dual
// Ey is primal dual primal
// the operation is thus dpd to pdp
auto constexpr p2dShift = primalToDual();
auto constexpr d2pShift = dualToPrimal();

if constexpr (dimension == 1)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0}, 0.5};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift}, 0.5};
return std::array<WeightPoint<dimension>, 2>{P1, P2};
}
else if constexpr (dimension == 2)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0}, 0.25};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift, 0}, 0.25};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, p2dShift}, 0.25};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, p2dShift},
0.25};
return std::array<WeightPoint<dimension>, 4>{P1, P2, P3, P4};
}
else if constexpr (dimension == 3)
{
constexpr WeightPoint<dimension> P1{Point<int, dimension>{0, 0, 0}, 0.25};
constexpr WeightPoint<dimension> P2{Point<int, dimension>{d2pShift, 0, 0}, 0.25};
constexpr WeightPoint<dimension> P3{Point<int, dimension>{0, p2dShift, 0}, 0.25};
constexpr WeightPoint<dimension> P4{Point<int, dimension>{d2pShift, p2dShift, 0},
0.25};
constexpr WeightPoint<dimension> P5{Point<int, dimension>{0, 0, d2pShift}, 0.25};
constexpr WeightPoint<dimension> P6{Point<int, dimension>{d2pShift, 0, d2pShift},
0.25};
constexpr WeightPoint<dimension> P7{Point<int, dimension>{0, p2dShift, d2pShift},
0.25};
constexpr WeightPoint<dimension> P8{
Point<int, dimension>{d2pShift, p2dShift, d2pShift}, 0.25};
return std::array<WeightPoint<dimension>, 8>{P1, P2, P3, P4, P5, P6, P7, P8};
}
}


NO_DISCARD auto static constexpr BzToEy()
Expand Down
13 changes: 9 additions & 4 deletions src/core/data/vecfield/vecfield.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,16 @@ namespace core
Vavg.getComponent(Component::Z));
}

// template<std::size_t... Index>
// NO_DISCARD auto norm(std::index_sequence<Index...>) const
// template<typename VecField, typename... Index>
// NO_DISCARD auto norm(VecField const& vf, Index... index)
// {
// std::sqrt(std::accumulate(std::begin(components_), std::end(components_), Type{0},
// [](auto acc, auto const& c) { return acc + c * c; }));
// using Type = typename VecField::value_type;
// std::sqrt(std::accumulate(std::begin(vf.components()), std::end(vf.components()),
// Type{0},
// [&](auto acc, auto const& c) {
// auto const v = c(index...);
// return acc(index...) + v * v;
// }));
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

struct VecFieldNames
Expand Down
43 changes: 39 additions & 4 deletions src/core/numerics/ohm/ohm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

#include "initializer/data_provider.hpp"

#include <numeric>


namespace PHARE::core
{
Expand Down Expand Up @@ -52,10 +54,11 @@ class Ohm : public LayoutHolder<GridLayout>




private:
double const eta_;
double const nu_;

private:
template<typename VecField, typename Field>
struct OhmPack
{
Expand Down Expand Up @@ -281,6 +284,7 @@ class Ohm : public LayoutHolder<GridLayout>
return -gradPOnEy / nOnEy;
}
else
#include <cmath>
{
return 0.;
}
Expand Down Expand Up @@ -334,11 +338,42 @@ class Ohm : public LayoutHolder<GridLayout>



template<auto component, typename VecField>
auto hyperresistive_(VecField const& J, VecField const& B, VecField const& n,
template<auto component, typename VecField, typename Field>
auto hyperresistive_(VecField const& J, VecField const& B, Field const& n,
MeshIndex<VecField::dimension> index) const
{ // TODO : https://github.com/PHAREHUB/PHARE/issues/3
return -nu_ * layout_->laplacian(J(component), index);

double const dl2{std::accumulate(std::begin(layout_->meshSize()),
std::end(layout_->meshSize()), 0.,
[](double acc, double d) { return acc + d * d; })};

if constexpr (component == Component::X)
{
auto const BxOnE = GridLayout::project(B(Component::X), index, GridLayout::BxToEx());
auto const ByOnE = GridLayout::project(B(Component::Y), index, GridLayout::ByToEx());
auto const BzOnE = GridLayout::project(B(Component::Z), index, GridLayout::BzToEx());
auto const nOnE = GridLayout::project(n, index, GridLayout::momentsToEx());
auto b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE);
return -nu_ * b / nOnE * dl2 * layout_->laplacian(J(component), index);
}
if constexpr (component == Component::Y)
{
auto const BxOnE = GridLayout::project(B(Component::X), index, GridLayout::BxToEy());
auto const ByOnE = GridLayout::project(B(Component::Y), index, GridLayout::ByToEy());
auto const BzOnE = GridLayout::project(B(Component::Z), index, GridLayout::BzToEy());
auto const nOnE = GridLayout::project(n, index, GridLayout::momentsToEy());
auto b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE);
return -nu_ * b / nOnE * dl2 * layout_->laplacian(J(component), index);
}
if constexpr (component == Component::Z)
{
auto const BxOnE = GridLayout::project(B(Component::X), index, GridLayout::BxToEz());
auto const ByOnE = GridLayout::project(B(Component::Y), index, GridLayout::ByToEz());
auto const BzOnE = GridLayout::project(B(Component::Z), index, GridLayout::BzToEz());
auto const nOnE = GridLayout::project(n, index, GridLayout::momentsToEz());
auto b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE);
return -nu_ * b / nOnE * dl2 * layout_->laplacian(J(component), index);
}
}
};

Expand Down

0 comments on commit 9a12fd0

Please sign in to comment.