diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 4d4df6deb..6b7356b28 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -236,8 +236,10 @@ def as_paths(rb): ) # integrator.h might want some looking at add_string("simulation/algo/ion_updater/pusher/name", simulation.particle_pusher) + add_double("simulation/algo/ohm/resistivity", simulation.resistivity) add_double("simulation/algo/ohm/hyper_resistivity", simulation.hyper_resistivity) + add_string("simulation/algo/ohm/hyper_mode", simulation.hyper_mode) # load balancer block start lb = simulation.load_balancer or LoadBalancer(active=False, _register=False) diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 5c36f67e3..5fb41ecc3 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -604,7 +604,11 @@ def check_hyper_resistivity(**kwargs): if hyper_resistivity < 0.0: raise ValueError("Error: hyper_resistivity should not be negative") - return hyper_resistivity + hyper_mode = kwargs.get("hyper_mode", "constant") + if hyper_mode not in ["constant", "spatial"]: + raise ValueError("Error: hyper_mode should be 'constant' or 'spatial'") + + return hyper_resistivity, hyper_mode def check_clustering(**kwargs): @@ -647,6 +651,7 @@ def wrapper(simulation_object, **kwargs): "diag_options", "resistivity", "hyper_resistivity", + "hyper_mode", "strict", "restart_options", "tag_buffer", @@ -718,7 +723,9 @@ def wrapper(simulation_object, **kwargs): kwargs["resistivity"] = check_resistivity(**kwargs) - kwargs["hyper_resistivity"] = check_hyper_resistivity(**kwargs) + nu, hyper_mode = check_hyper_resistivity(**kwargs) + kwargs["hyper_resistivity"] = nu + kwargs["hyper_mode"] = hyper_mode kwargs["dry_run"] = kwargs.get( "dry_run", os.environ.get("PHARE_DRY_RUN", "0") == "1" diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp index 4faf3bb39..0efd3f795 100644 --- a/src/amr/resources_manager/amr_utils.hpp +++ b/src/amr/resources_manager/amr_utils.hpp @@ -182,7 +182,8 @@ namespace amr nbrCell[iDim] = static_cast(domain.numberCells(iDim)); } - return GridLayoutT{dl, nbrCell, origin, amr::Box{domain}}; + auto lvlNbr = patch.getPatchLevelNumber(); + return GridLayoutT{dl, nbrCell, origin, amr::Box{domain}, lvlNbr}; } inline auto to_string(SAMRAI::hier::GlobalId const& id) diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 59bd1ef8b..74297a864 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -112,7 +112,7 @@ namespace core GridLayout(std::array const& meshSize, std::array const& nbrCells, Point const& origin, - Box AMRBox = Box{}) + Box AMRBox = Box{}, int level_number = 0) : meshSize_{meshSize} , origin_{origin} , nbrPhysicalCells_{nbrCells} @@ -120,6 +120,7 @@ namespace core , physicalEndIndexTable_{initPhysicalEnd_()} , ghostEndIndexTable_{initGhostEnd_()} , AMRBox_{AMRBox} + , levelNumber_{level_number} { if (AMRBox_.isEmpty()) { @@ -1073,6 +1074,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 +1090,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 +1098,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 +1122,8 @@ namespace core */ NO_DISCARD auto static constexpr ByToEz() { return GridLayoutImpl::ByToEz(); } + NO_DISCARD auto static constexpr BzToEz() { return GridLayoutImpl::BzToEz(); } + /** @@ -1165,6 +1171,7 @@ namespace core evalOnBox_(field, fn, indices); } + auto levelNumber() const { return levelNumber_; } private: template @@ -1508,6 +1515,8 @@ namespace core // arrays will be accessed with [primal] and [dual] indexes. constexpr static std::array nextIndexTable_{{nextPrimal_(), nextDual_()}}; constexpr static std::array prevIndexTable_{{prevPrimal_(), prevDual_()}}; + + int levelNumber_ = 0; }; 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/tensorfield/tensorfield.hpp b/src/core/data/tensorfield/tensorfield.hpp index 991e18d42..ffc6bed92 100644 --- a/src/core/data/tensorfield/tensorfield.hpp +++ b/src/core/data/tensorfield/tensorfield.hpp @@ -82,15 +82,13 @@ class TensorField NO_DISCARD auto getCompileTimeResourcesViewList() { - return for_N([&](auto i) -> auto& { - return components_[i]; - }); + return for_N( + [&](auto i) -> auto& { return components_[i]; }); } NO_DISCARD auto getCompileTimeResourcesViewList() const { - return for_N([&](auto i) -> auto& { - return components_[i]; - }); + return for_N( + [&](auto i) -> auto& { return components_[i]; }); } diff --git a/src/core/data/vecfield/vecfield.hpp b/src/core/data/vecfield/vecfield.hpp index 1fa9f4fec..746bac66d 100644 --- a/src/core/data/vecfield/vecfield.hpp +++ b/src/core/data/vecfield/vecfield.hpp @@ -32,6 +32,7 @@ namespace core } + struct VecFieldNames { std::string vecName; diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index f0e5fefa5..070364320 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -1,19 +1,19 @@ #ifndef PHARE_OHM_HPP #define PHARE_OHM_HPP -#include -#include #include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/grid/gridlayout.hpp" #include "core/data/grid/gridlayout_utils.hpp" #include "core/data/vecfield/vecfield_component.hpp" #include "initializer/data_provider.hpp" + namespace PHARE::core { +enum class HyperMode { constant, spatial }; + template class Ohm : public LayoutHolder { @@ -24,6 +24,9 @@ class Ohm : public LayoutHolder explicit Ohm(PHARE::initializer::PHAREDict const& dict) : 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} { } @@ -52,10 +55,12 @@ class Ohm : public LayoutHolder + +private: double const eta_; double const nu_; + HyperMode const hyper_mode; -private: template struct OhmPack { @@ -76,7 +81,7 @@ class Ohm : public LayoutHolder Exyz(ijk...) = ideal_(Ve, B, {ijk...}) // + pressure_(n, Pe, {ijk...}) // + resistive_(J, {ijk...}) // - + hyperresistive_(J, {ijk...}); + + hyperresistive_(J, B, n, {ijk...}); } @@ -331,14 +336,57 @@ class Ohm : public LayoutHolder } } - + template + auto hyperresistive_(VecField const& J, VecField const& B, Field const& n, + MeshIndex index) const + { + if (hyper_mode == HyperMode::constant) + return constant_hyperresistive_(J, index); + else if (hyper_mode == HyperMode::spatial) + return spatial_hyperresistive_(J, B, n, index); + else // should not happen but otherwise -Wreturn-type fails with Werror + throw std::runtime_error("Error - Ohm - unknown hyper_mode"); + } template - auto hyperresistive_(VecField const& J, MeshIndex index) const + auto constant_hyperresistive_(VecField const& J, MeshIndex index) const { // TODO : https://github.com/PHAREHUB/PHARE/issues/3 return -nu_ * layout_->laplacian(J(component), index); } + + + template + 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 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 * lvlCoeff * layout_->laplacian(J(component), index); + }; + + if constexpr (component == Component::X) + { + return computeHR(GridLayout::BxToEx(), GridLayout::ByToEx(), GridLayout::BzToEx(), + GridLayout::momentsToEx()); + } + if constexpr (component == Component::Y) + { + return computeHR(GridLayout::BxToEy(), GridLayout::ByToEy(), GridLayout::BzToEy(), + GridLayout::momentsToEy()); + } + if constexpr (component == Component::Z) + { + return computeHR(GridLayout::BxToEz(), GridLayout::ByToEz(), GridLayout::BzToEz(), + GridLayout::momentsToEz()); + } + } };