hyperresistivity can depend on level and local B/n

PHAREHUB · Oct 12, 2024 · 31cdd1b · 31cdd1b
1 parent ecef2fe
commit 31cdd1b
Show file tree

Hide file tree

Showing 8 changed files with 199 additions and 15 deletions.
diff --git 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
@@ -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
@@ -182,7 +182,8 @@ namespace amr
             nbrCell[iDim] = static_cast<std::uint32_t>(domain.numberCells(iDim));
         }
 
-        return GridLayoutT{dl, nbrCell, origin, amr::Box<int, dimension>{domain}};
+        auto lvlNbr = patch.getPatchLevelNumber();
+        return GridLayoutT{dl, nbrCell, origin, amr::Box<int, dimension>{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
@@ -112,14 +112,15 @@ namespace core
         GridLayout(std::array<double, dimension> const& meshSize,
                    std::array<std::uint32_t, dimension> const& nbrCells,
                    Point<double, dimension> const& origin,
-                   Box<int, dimension> AMRBox = Box<int, dimension>{})
+                   Box<int, dimension> AMRBox = Box<int, dimension>{}, int level_number = 0)
             : meshSize_{meshSize}
             , origin_{origin}
             , nbrPhysicalCells_{nbrCells}
             , physicalStartIndexTable_{initPhysicalStart_()}
             , 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,13 +1090,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
@@ -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<typename Field, typename IndicesFn, typename Fn>
@@ -1508,6 +1515,8 @@ namespace core
         // arrays will be accessed with [primal] and [dual] indexes.
         constexpr static std::array<int, 2> nextIndexTable_{{nextPrimal_(), nextDual_()}};
         constexpr static std::array<int, 2> prevIndexTable_{{prevPrimal_(), prevDual_()}};
+
+        int levelNumber_ = 0;
     };
 
 

diff --git 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<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()
@@ -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()
@@ -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()

diff --git 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<N, for_N_R_mode::forward_tuple>([&](auto i) -> auto& {
-            return components_[i];
-        });
+        return for_N<N, for_N_R_mode::forward_tuple>(
+            [&](auto i) -> auto& { return components_[i]; });
     }
     NO_DISCARD auto getCompileTimeResourcesViewList() const
     {
-        return for_N<N, for_N_R_mode::forward_tuple>([&](auto i) -> auto& {
-            return components_[i];
-        });
+        return for_N<N, for_N_R_mode::forward_tuple>(
+            [&](auto i) -> auto& { return components_[i]; });
     }
 
 

diff --git 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
@@ -1,19 +1,19 @@
 #ifndef PHARE_OHM_HPP
 #define PHARE_OHM_HPP
 
-#include <cstddef>
-#include <iostream>
 
 #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<typename GridLayout>
 class Ohm : public LayoutHolder<GridLayout>
 {
@@ -24,6 +24,9 @@ class Ohm : public LayoutHolder<GridLayout>
     explicit Ohm(PHARE::initializer::PHAREDict const& dict)
         : eta_{dict["resistivity"].template to<double>()}
         , nu_{dict["hyper_resistivity"].template to<double>()}
+        , hyper_mode{cppdict::get_value(dict, "hyper_mode", std::string{"constant"}) == "constant"
+                         ? HyperMode::constant
+                         : HyperMode::spatial}
     {
     }
 
@@ -52,10 +55,12 @@ class Ohm : public LayoutHolder<GridLayout>
 
 
 
+
+private:
     double const eta_;
     double const nu_;
+    HyperMode const hyper_mode;
 
-private:
     template<typename VecField, typename Field>
     struct OhmPack
     {
@@ -76,7 +81,7 @@ class Ohm : public LayoutHolder<GridLayout>
         Exyz(ijk...) = ideal_<Tag>(Ve, B, {ijk...})      //
                        + pressure_<Tag>(n, Pe, {ijk...}) //
                        + resistive_<Tag>(J, {ijk...})    //
-                       + hyperresistive_<Tag>(J, {ijk...});
+                       + hyperresistive_<Tag>(J, B, n, {ijk...});
     }
 
 
@@ -281,6 +286,7 @@ class Ohm : public LayoutHolder<GridLayout>
                 return -gradPOnEy / nOnEy;
             }
             else
+#include <cmath>
             {
                 return 0.;
             }
@@ -331,14 +337,57 @@ class Ohm : public LayoutHolder<GridLayout>
         }
     }
 
-
+    template<auto component, typename VecField, typename Field>
+    auto hyperresistive_(VecField const& J, VecField const& B, Field const& n,
+                         MeshIndex<VecField::dimension> index) const
+    {
+        if (hyper_mode == HyperMode::constant)
+            return constant_hyperresistive_<component>(J, index);
+        else if (hyper_mode == HyperMode::spatial)
+            return spatial_hyperresistive_<component>(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 component, typename VecField>
-    auto hyperresistive_(VecField const& J, MeshIndex<VecField::dimension> index) const
+    auto constant_hyperresistive_(VecField const& J, MeshIndex<VecField::dimension> index) const
     { // TODO : https://github.com/PHAREHUB/PHARE/issues/3
         return -nu_ * layout_->laplacian(J(component), index);
     }
+
+
+    template<auto component, typename VecField, typename Field>
+    auto spatial_hyperresistive_(VecField const& J, VecField const& B, Field const& n,
+                                 MeshIndex<VecField::dimension> 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());
+        }
+    }
 };