Skip to content

Commit

Permalink
spatially variable hyperresistivity (#893)
Browse files Browse the repository at this point in the history 
* silent werror on ubuntu builds

* debug03samrai

* Update cmake_ubuntu.yml

* hyperresistivity can depend on level and local B/n

* hypermode check

* less typo

---------

Co-authored-by: PhilipDeegan <[email protected]>
  • Loading branch information
@nicolasaunai @PhilipDeegan
nicolasaunai and PhilipDeegan authored Oct 18, 2024
1 parent a69bcf0 commit af3f2be
Show file tree
Hide file tree
Showing 8 changed files with 205 additions and 17 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
11 changes: 9 additions & 2 deletions pyphare/pyphare/pharein/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -647,6 +651,7 @@ def wrapper(simulation_object, **kwargs):
"diag_options",
"resistivity",
"hyper_resistivity",
"hyper_mode",
"strict",
"restart_options",
"tag_buffer",
Expand Down Expand Up @@ -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"
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
1 change: 1 addition & 0 deletions src/core/data/vecfield/vecfield.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ namespace core
}



struct VecFieldNames
{
std::string vecName;
Expand Down
Loading

0 comments on commit af3f2be

Please sign in to comment.