Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

absolute tagging #902

Closed
wants to merge 11 commits into from
Closed

absolute tagging #902

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
6822967
wip
nicolasaunai Sep 16, 2024
0f755f9
variable nu
nicolasaunai Sep 18, 2024
1fa03ba
add hyper_mode
nicolasaunai Sep 20, 2024
85df275
Wreturn-type
nicolasaunai Sep 20, 2024
f5f2f05
addstring
nicolasaunai Sep 20, 2024
e513179
refacrabbit
nicolasaunai Sep 20, 2024
d03d643
level coef in hypermode
nicolasaunai Sep 27, 2024
08ce776
accept hyper_mode kw
nicolasaunai Oct 1, 2024
b7b6db9
level coef in right order
nicolasaunai Oct 4, 2024
4c99903
keep particle overallocation on restore state (#897)
PhilipDeegan Oct 3, 2024
fef5511
tagging better
nicolasaunai Oct 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions pyphare/pyphare/pharein/__init__.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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};
Comment on lines +185 to +186
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🛠️ Refactor suggestion

Consider improving backwards compatibility and documentation

The addition of the patch level number (lvlNbr) to the GridLayoutT constructor is a valuable change that can enhance level-specific operations and debugging. However, this modification has some implications:

  1. The function signature change may break existing code that uses layoutFromPatch.
  2. The function's documentation hasn't been updated to reflect this new parameter.

To address these issues:

  1. Consider using a default parameter to maintain backwards compatibility:

    NO_DISCARD GridLayoutT layoutFromPatch(SAMRAI::hier::Patch const& patch, int lvlNbr = -1)

    If lvlNbr is not provided, you can still obtain it from patch.getPatchLevelNumber().

  2. Update the function documentation to describe the new lvlNbr parameter and its purpose.

Here's a suggested implementation with improved documentation:

/**
 * @brief Create a GridLayoutT from a SAMRAI::hier::Patch
 * @param patch The SAMRAI patch to create the layout from
 * @param lvlNbr The patch level number (optional, default: -1)
 * @return A GridLayoutT instance representing the patch
 */
template<typename GridLayoutT>
NO_DISCARD GridLayoutT layoutFromPatch(SAMRAI::hier::Patch const& patch, int lvlNbr = -1)
{
    // ... (existing code)

    if (lvlNbr == -1)
    {
        lvlNbr = patch.getPatchLevelNumber();
    }

    return GridLayoutT{dl, nbrCell, origin, amr::Box<int, dimension>{domain}, lvlNbr};
}

This approach maintains backwards compatibility while allowing explicit specification of the level number when needed.

}

inline auto to_string(SAMRAI::hier::GlobalId const& id)
Expand Down
28 changes: 19 additions & 9 deletions src/amr/solvers/solver_ppc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,22 @@ class SolverPPC : public ISolver<AMR_Types>
std::unordered_map<std::string, ParticleArray> tmpDomain;
std::unordered_map<std::string, ParticleArray> patchGhost;

template<typename Map>
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

Expand Down Expand Up @@ -206,15 +222,9 @@ void SolverPPC<HybridModel, AMR_Types>::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());
}
}
}
Expand Down
7 changes: 4 additions & 3 deletions src/amr/tagging/default_hybrid_tagger_strategy.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,9 +108,10 @@ void DefaultHybridTaggerStrategy<HybridModel>::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));
Comment on lines +111 to +114
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🛠️ Refactor suggestion

Approved: Changes align with PR objectives and improve formula consistency.

The modifications to the field_diff lambda function address the PR objective of incorporating absolute values in the tagging formula. This should help in refining the two current layers uniformly across different thresholds.

A minor optimization suggestion:

Consider factoring out the common std::abs() call in the denominator of the second tuple element:

 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));
+                            / (1 + std::abs(F(ix, iy + 1) - F(ix, iy))));

This change makes the second tuple element consistent with the first one and potentially improves readability.

📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
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));
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))
/ (1 + std::abs(F(ix, iy + 1) - F(ix, iy))));

};

auto const& [Bx_x, Bx_y] = field_diff(Bx);
Expand Down
11 changes: 10 additions & 1 deletion src/core/data/grid/gridlayout.hpp
View file
Open in desktop
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)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue

Ensure consistent parameter type for level_number

The new constructor parameter int level_number = 0 introduces the level_number as an int. Since level_number represents a non-negative value (level index), consider using std::size_t or unsigned int for consistency and to prevent negative values.

Apply this diff to change the parameter type:

-Box<int, dimension> AMRBox = Box<int, dimension>{}, int level_number = 0)
+Box<int, dimension> AMRBox = Box<int, dimension>{}, std::size_t level_number = 0)
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
Box<int, dimension> AMRBox = Box<int, dimension>{}, int level_number = 0)
Box<int, dimension> AMRBox = Box<int, dimension>{}, std::size_t 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;
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🛠️ Refactor suggestion

Consider using std::size_t for levelNumber_

Since levelNumber_ represents a non-negative level index, using an unsigned type like std::size_t can prevent negative values and improve type safety.

Apply this diff:

-int levelNumber_;
+std::size_t levelNumber_;
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
int levelNumber_ = 0;
std::size_t levelNumber_ = 0;

};


Expand Down
122 changes: 122 additions & 0 deletions src/core/data/grid/gridlayoutimplyee.hpp
View file
Open in desktop
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};
}
}
Comment on lines +688 to +728
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🛠️ Refactor suggestion

Consider refactoring BxToEx() to reduce code duplication

The BxToEx() method contains similar code blocks for different dimensions (dimension == 1, dimension == 2, and dimension == 3). This can lead to code duplication and increase maintenance effort. To enhance maintainability and readability, consider refactoring by creating a helper function or using templates to handle the common logic across dimensions.



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
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,17 @@
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;
// }));
// }
Comment on lines +34 to +44

Check notice

Code scanning / CodeQL

Commented-out code Note

This comment appears to contain commented-out code.

struct VecFieldNames
{
Expand Down
Loading
Loading