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..1f4fea60d 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -647,6 +647,7 @@ def wrapper(simulation_object, **kwargs): "diag_options", "resistivity", "hyper_resistivity", + "hyper_mode", "strict", "restart_options", "tag_buffer", @@ -719,6 +720,7 @@ def wrapper(simulation_object, **kwargs): kwargs["resistivity"] = check_resistivity(**kwargs) kwargs["hyper_resistivity"] = check_hyper_resistivity(**kwargs) + kwargs["hyper_mode"] = kwargs.get("hyper_mode", "constant") 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/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index eb1f79fe8..3eb79f496 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -152,6 +152,22 @@ class SolverPPC : public ISolver std::unordered_map tmpDomain; std::unordered_map patchGhost; + template + static void add_to(Map& map, std::string const& key, ParticleArray const& ps) + { + // vector copy drops the capacity (over allocation of the source) + // we want to keep the overallocation somewhat - how much to be assessed + ParticleArray empty{ps.box()}; + + if (!map.count(key)) + map.emplace(key, empty); + else + map.at(key) = empty; + + auto& v = map.at(key); + v.reserve(ps.capacity()); + v.replace_from(ps); + } }; // end solverPPC @@ -206,15 +222,9 @@ void SolverPPC::saveState_(level_t& level, ModelViews_t& ss << state.patch->getGlobalId(); for (auto& pop : state.ions) { - auto key = ss.str() + "_" + pop.name(); - if (!tmpDomain.count(key)) - tmpDomain.emplace(key, pop.domainParticles()); - else - tmpDomain.at(key) = pop.domainParticles(); - if (!patchGhost.count(key)) - patchGhost.emplace(key, pop.patchGhostParticles()); - else - patchGhost.at(key) = pop.patchGhostParticles(); + std::string const key = ss.str() + "_" + pop.name(); + add_to(tmpDomain, key, pop.domainParticles()); + add_to(patchGhost, key, pop.patchGhostParticles()); } } } diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp index 2fb3ea9df..ffac48b93 100644 --- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp +++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp @@ -108,9 +108,10 @@ void DefaultHybridTaggerStrategy::tag(HybridModel& model, { auto field_diff = [&](auto const& F) // { - return std::make_tuple( - std::abs((F(ix + 2, iy) - F(ix, iy)) / (1 + F(ix + 1, iy) - F(ix, iy))), - std::abs((F(ix, iy + 2) - F(ix, iy)) / (F(ix, iy + 1) - F(ix, iy) + 1))); + return std::make_tuple(std::abs((F(ix + 2, iy) - F(ix, iy)) + / (1 + std::abs(F(ix + 1, iy) - F(ix, iy)))), + std::abs(F(ix, iy + 2) - F(ix, iy)) + / (std::abs(F(ix, iy + 1) - F(ix, iy)) + 1)); }; auto const& [Bx_x, Bx_y] = field_diff(Bx); 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..ad0d2357b 100644 --- a/src/core/data/vecfield/vecfield.hpp +++ b/src/core/data/vecfield/vecfield.hpp @@ -31,6 +31,17 @@ namespace core Vavg.getComponent(Component::Z)); } + // template + // NO_DISCARD auto norm(VecField const& vf, Index... index) + // { + // 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 f0e5fefa5..b4bc34ddf 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -11,9 +11,13 @@ #include "initializer/data_provider.hpp" +#include + namespace PHARE::core { +enum class HyperMode { constant, spatial }; + template class Ohm : public LayoutHolder { @@ -24,6 +28,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 +59,12 @@ class Ohm : public LayoutHolder + +private: double const eta_; double const nu_; + HyperMode const hyper_mode; -private: template struct OhmPack { @@ -76,7 +85,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...}); } @@ -281,6 +290,7 @@ class Ohm : public LayoutHolder return -gradPOnEy / nOnEy; } else +#include { return 0.; } @@ -331,14 +341,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 + return 0.; + } 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()); + } + } };