From 6822967e5b72fe18e5d874d2f8c9253bbbe99d69 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Mon, 16 Sep 2024 10:29:58 +0200 Subject: [PATCH 01/11] wip --- src/core/data/tensorfield/tensorfield.hpp | 10 ++++------ src/core/data/vecfield/vecfield.hpp | 6 ++++++ src/core/numerics/ohm/ohm.hpp | 5 +++-- 3 files changed, 13 insertions(+), 8 deletions(-) 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..a0164ba42 100644 --- a/src/core/data/vecfield/vecfield.hpp +++ b/src/core/data/vecfield/vecfield.hpp @@ -31,6 +31,12 @@ namespace core Vavg.getComponent(Component::Z)); } + // template + // NO_DISCARD auto norm(std::index_sequence) const + // { + // std::sqrt(std::accumulate(std::begin(components_), std::end(components_), Type{0}, + // [](auto acc, auto const& c) { return acc + c * c; })); + // } struct VecFieldNames { diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index f0e5fefa5..5abc26501 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -76,7 +76,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...}); } @@ -335,7 +335,8 @@ class Ohm : public LayoutHolder template - auto hyperresistive_(VecField const& J, MeshIndex index) const + auto hyperresistive_(VecField const& J, VecField const& B, VecField const& n, + MeshIndex index) const { // TODO : https://github.com/PHAREHUB/PHARE/issues/3 return -nu_ * layout_->laplacian(J(component), index); } From 0f755f92edd24f5f8500ab032a39704784c8e839 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Wed, 18 Sep 2024 16:49:33 +0200 Subject: [PATCH 02/11] 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); + } } }; From 1fa03ba8355fbaed174adf3b8b3276f74d1f8ddd Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 20 Sep 2024 09:32:04 +0200 Subject: [PATCH 03/11] add hyper_mode --- pyphare/pyphare/pharein/__init__.py | 2 ++ pyphare/pyphare/pharein/simulation.py | 1 + src/core/numerics/ohm/ohm.hpp | 26 +++++++++++++++++++++++--- 3 files changed, 26 insertions(+), 3 deletions(-) diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 4d4df6deb..707edb31c 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_double("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..022410f1c 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -719,6 +719,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/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index e933e7b33..5b777ff11 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -16,6 +16,8 @@ namespace PHARE::core { +enum class HyperMode { constant, spatial }; + template class Ohm : public LayoutHolder { @@ -26,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} { } @@ -58,6 +63,7 @@ class Ohm : public LayoutHolder private: double const eta_; double const nu_; + HyperMode const hyper_mode; template struct OhmPack @@ -335,14 +341,28 @@ 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); + } + + + template + 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 double const dl2{std::accumulate(std::begin(layout_->meshSize()), std::end(layout_->meshSize()), 0., [](double acc, double d) { return acc + d * d; })}; From 85df275839c509180415e54f0a18e28235a1af6b Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 20 Sep 2024 09:40:43 +0200 Subject: [PATCH 04/11] Wreturn-type --- src/core/numerics/ohm/ohm.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 5b777ff11..eb3f7377e 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -349,6 +349,8 @@ class Ohm : public LayoutHolder 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.; } From f5f2f05f07f5decd159d0f1b4b84c11b8d784ec6 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 20 Sep 2024 10:10:47 +0200 Subject: [PATCH 05/11] addstring --- pyphare/pyphare/pharein/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 707edb31c..6b7356b28 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -239,7 +239,7 @@ def as_paths(rb): add_double("simulation/algo/ohm/resistivity", simulation.resistivity) add_double("simulation/algo/ohm/hyper_resistivity", simulation.hyper_resistivity) - add_double("simulation/algo/ohm/hyper_mode", simulation.hyper_mode) + add_string("simulation/algo/ohm/hyper_mode", simulation.hyper_mode) # load balancer block start lb = simulation.load_balancer or LoadBalancer(active=False, _register=False) From e5131794f53b219e921513cf0a23ad54c46ab0ea Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 20 Sep 2024 17:05:32 +0200 Subject: [PATCH 06/11] refacrabbit --- src/core/numerics/ohm/ohm.hpp | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index eb3f7377e..8d3c60927 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -369,32 +369,29 @@ class Ohm : public LayoutHolder 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 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 * dl2 * 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) { - 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); + return computeHR(GridLayout::BxToEy(), GridLayout::ByToEy(), GridLayout::BzToEy(), + GridLayout::momentsToEy()); } 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); + return computeHR(GridLayout::BxToEz(), GridLayout::ByToEz(), GridLayout::BzToEz(), + GridLayout::momentsToEz()); } } }; From d03d643dcd536947bed77aaf257415d93000d870 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 27 Sep 2024 15:42:37 +0200 Subject: [PATCH 07/11] level coef in hypermode --- src/amr/resources_manager/amr_utils.hpp | 3 ++- src/core/data/grid/gridlayout.hpp | 6 +++++- src/core/numerics/ohm/ohm.hpp | 6 ++---- 3 files changed, 9 insertions(+), 6 deletions(-) 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 9ab41c674..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()) { @@ -1170,6 +1171,7 @@ namespace core evalOnBox_(field, fn, indices); } + auto levelNumber() const { return levelNumber_; } private: template @@ -1513,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/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 8d3c60927..39a27172e 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -365,9 +365,7 @@ class Ohm : public LayoutHolder auto spatial_hyperresistive_(VecField const& J, VecField const& B, Field const& n, MeshIndex index) const { // TODO : https://github.com/PHAREHUB/PHARE/issues/3 - double const dl2{std::accumulate(std::begin(layout_->meshSize()), - std::end(layout_->meshSize()), 0., - [](double acc, double d) { return acc + d * d; })}; + auto const lvlCoeff = 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); @@ -375,7 +373,7 @@ class Ohm : public LayoutHolder 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 * dl2 * layout_->laplacian(J(component), index); + return -nu_ * b / nOnE * lvlCoeff * layout_->laplacian(J(component), index); }; if constexpr (component == Component::X) From 08ce776b1fd24e65a624e2efe473f648f0b624fb Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 1 Oct 2024 10:40:08 +0200 Subject: [PATCH 08/11] accept hyper_mode kw --- pyphare/pyphare/pharein/simulation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 022410f1c..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", From b7b6db9356bd1b5e527de1637ccb6c75810b5382 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 4 Oct 2024 16:47:46 +0200 Subject: [PATCH 09/11] level coef in right order --- src/core/numerics/ohm/ohm.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 39a27172e..b4bc34ddf 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -365,7 +365,7 @@ class Ohm : public LayoutHolder 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 = std::pow(4, layout_->levelNumber()); + 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); From 4c9990374ad718311b3879e14d319363d0b6b05c Mon Sep 17 00:00:00 2001 From: PhilipDeegan Date: Thu, 3 Oct 2024 17:58:59 +0200 Subject: [PATCH 10/11] keep particle overallocation on restore state (#897) --- src/amr/solvers/solver_ppc.hpp | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) 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()); } } } From fef5511ac013872567f54e8eb5ee3af6b0f1038a Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 11 Oct 2024 09:35:04 +0200 Subject: [PATCH 11/11] tagging better --- src/amr/tagging/default_hybrid_tagger_strategy.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) 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);