From 9a12fd058064a6db55218e6735e8a46ed64acf88 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Wed, 18 Sep 2024 16:49:33 +0200 Subject: [PATCH] variable nu --- src/core/data/grid/gridlayout.hpp | 5 + src/core/data/grid/gridlayoutimplyee.hpp | 122 +++++++++++++++++++++++ src/core/data/vecfield/vecfield.hpp | 13 ++- src/core/numerics/ohm/ohm.hpp | 43 +++++++- 4 files changed, 175 insertions(+), 8 deletions(-) diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 59bd1ef8b..9ab41c674 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -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 @@ -1088,6 +1089,7 @@ 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. @@ -1095,6 +1097,7 @@ namespace core 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 @@ -1118,6 +1121,8 @@ namespace core */ NO_DISCARD auto static constexpr ByToEz() { return GridLayoutImpl::ByToEz(); } + NO_DISCARD auto static constexpr BzToEz() { return GridLayoutImpl::BzToEz(); } + /** diff --git a/src/core/data/grid/gridlayoutimplyee.hpp b/src/core/data/grid/gridlayoutimplyee.hpp index 031cf07b0..76926ce93 100644 --- a/src/core/data/grid/gridlayoutimplyee.hpp +++ b/src/core/data/grid/gridlayoutimplyee.hpp @@ -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 P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{p2dShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.25}; + constexpr WeightPoint P2{Point{0, d2pShift}, 0.25}; + constexpr WeightPoint P3{Point{p2dShift, 0}, 0.25}; + constexpr WeightPoint P4{Point{p2dShift, d2pShift}, + 0.25}; + return std::array, 4>{P1, P2, P3, P4}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.125}; + constexpr WeightPoint P2{Point{0, d2pShift, 0}, 0.125}; + constexpr WeightPoint P3{Point{p2dShift, 0, 0}, 0.125}; + constexpr WeightPoint P4{Point{p2dShift, d2pShift, 0}, + 0.125}; + + constexpr WeightPoint P5{Point{0, 0, d2pShift}, 0.125}; + constexpr WeightPoint P6{Point{0, d2pShift, d2pShift}, + 0.125}; + constexpr WeightPoint P7{Point{p2dShift, 0, d2pShift}, + 0.125}; + constexpr WeightPoint P8{ + Point{p2dShift, d2pShift, d2pShift}, 0.125}; + return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + } + } NO_DISCARD auto static constexpr ByToEx() @@ -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 P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{d2pShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.25}; + constexpr WeightPoint P2{Point{d2pShift, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, d2pShift}, 0.25}; + constexpr WeightPoint P4{Point{d2pShift, d2pShift}, + 0.25}; + return std::array, 4>{P1, P2, P3, P4}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.25}; + constexpr WeightPoint P2{Point{d2pShift, 0, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, d2pShift, 0}, 0.25}; + constexpr WeightPoint P4{Point{d2pShift, d2pShift, 0}, + 0.25}; + constexpr WeightPoint P5{Point{0, 0, p2dShift}, 0.25}; + constexpr WeightPoint P6{Point{d2pShift, 0, p2dShift}, + 0.25}; + constexpr WeightPoint P7{Point{0, d2pShift, p2dShift}, + 0.25}; + constexpr WeightPoint P8{ + Point{d2pShift, d2pShift, p2dShift}, 0.25}; + return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + } + } NO_DISCARD auto static constexpr ByToEz() @@ -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 P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{d2pShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.25}; + constexpr WeightPoint P2{Point{d2pShift, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, p2dShift}, 0.25}; + constexpr WeightPoint P4{Point{d2pShift, p2dShift}, + 0.25}; + return std::array, 4>{P1, P2, P3, P4}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.25}; + constexpr WeightPoint P2{Point{d2pShift, 0, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, p2dShift, 0}, 0.25}; + constexpr WeightPoint P4{Point{d2pShift, p2dShift, 0}, + 0.25}; + constexpr WeightPoint P5{Point{0, 0, d2pShift}, 0.25}; + constexpr WeightPoint P6{Point{d2pShift, 0, d2pShift}, + 0.25}; + constexpr WeightPoint P7{Point{0, p2dShift, d2pShift}, + 0.25}; + constexpr WeightPoint P8{ + Point{d2pShift, p2dShift, d2pShift}, 0.25}; + return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + } + } NO_DISCARD auto static constexpr BzToEy() diff --git a/src/core/data/vecfield/vecfield.hpp b/src/core/data/vecfield/vecfield.hpp index a0164ba42..ad0d2357b 100644 --- a/src/core/data/vecfield/vecfield.hpp +++ b/src/core/data/vecfield/vecfield.hpp @@ -31,11 +31,16 @@ namespace core Vavg.getComponent(Component::Z)); } - // template - // NO_DISCARD auto norm(std::index_sequence) const + // template + // 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; + // })); // } struct VecFieldNames diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 5abc26501..e933e7b33 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -11,6 +11,8 @@ #include "initializer/data_provider.hpp" +#include + namespace PHARE::core { @@ -52,10 +54,11 @@ class Ohm : public LayoutHolder + +private: double const eta_; double const nu_; -private: template struct OhmPack { @@ -281,6 +284,7 @@ class Ohm : public LayoutHolder return -gradPOnEy / nOnEy; } else +#include { return 0.; } @@ -334,11 +338,42 @@ class Ohm : public LayoutHolder - template - auto hyperresistive_(VecField const& J, VecField const& B, VecField const& n, + template + auto hyperresistive_(VecField const& J, VecField const& B, Field const& n, MeshIndex 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); + } } };