Skip to content

Commit

Permalink
hyperresistivity can depend on level and local B/n
Browse files Browse the repository at this point in the history
  • Loading branch information
@nicolasaunai
nicolasaunai committed Oct 12, 2024
1 parent ecef2fe commit 8c22d42
Show file tree
Hide file tree
Showing 8 changed files with 210 additions and 12 deletions.
2 changes: 2 additions & 0 deletions pyphare/pyphare/pharein/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions pyphare/pyphare/pharein/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -647,6 +647,7 @@ def wrapper(simulation_object, **kwargs):
"diag_options",
"resistivity",
"hyper_resistivity",
"hyper_mode",
"strict",
"restart_options",
"tag_buffer",
Expand Down Expand Up @@ -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"
Expand Down
3 changes: 2 additions & 1 deletion src/amr/resources_manager/amr_utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
11 changes: 10 additions & 1 deletion src/core/data/grid/gridlayout.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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())
{
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -1118,6 +1122,8 @@ namespace core
*/
NO_DISCARD auto static constexpr ByToEz() { return GridLayoutImpl::ByToEz(); }

NO_DISCARD auto static constexpr BzToEz() { return GridLayoutImpl::BzToEz(); }



/**
Expand Down Expand Up @@ -1165,6 +1171,7 @@ namespace core
evalOnBox_(field, fn, indices);
}

auto levelNumber() const { return levelNumber_; }

private:
template<typename Field, typename IndicesFn, typename Fn>
Expand Down Expand Up @@ -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;
};


Expand Down
122 changes: 122 additions & 0 deletions src/core/data/grid/gridlayoutimplyee.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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()
Expand Down
10 changes: 4 additions & 6 deletions src/core/data/tensorfield/tensorfield.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]; });
}


Expand Down
11 changes: 11 additions & 0 deletions src/core/data/vecfield/vecfield.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,17 @@ namespace core
Vavg.getComponent(Component::Z));
}

// template<typename VecField, typename... Index>
// 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;
// }));
// }

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

struct VecFieldNames
{
Expand Down
61 changes: 57 additions & 4 deletions src/core/numerics/ohm/ohm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,13 @@

#include "initializer/data_provider.hpp"

#include <numeric>


namespace PHARE::core
{
enum class HyperMode { constant, spatial };

template<typename GridLayout>
class Ohm : public LayoutHolder<GridLayout>
{
Expand All @@ -24,6 +28,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}
{
}

Expand Down Expand Up @@ -52,10 +59,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
{
Expand All @@ -76,7 +85,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...});
}


Expand Down Expand Up @@ -281,6 +290,7 @@ class Ohm : public LayoutHolder<GridLayout>
return -gradPOnEy / nOnEy;
}
else
#include <cmath>
{
return 0.;
}
Expand Down Expand Up @@ -331,14 +341,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
return 0.;
}


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());
}
}
};


Expand Down

0 comments on commit 8c22d42

Please sign in to comment.