From b5faac985a9efb90592f265ca8823269bba6ce11 Mon Sep 17 00:00:00 2001 From: PhilipDeegan Date: Mon, 23 Mar 2026 22:37:03 +0100 Subject: [PATCH 1/3] drop ppc model view --- src/amr/multiphysics_integrator.hpp | 15 +- .../resources_manager/resources_manager.hpp | 24 +- src/amr/solvers/solver.hpp | 46 +-- src/amr/solvers/solver_field_evolvers.hpp | 117 +++++++ src/amr/solvers/solver_mhd.hpp | 34 +- src/amr/solvers/solver_ppc.hpp | 212 ++++++------- src/amr/solvers/solver_ppc_model_view.hpp | 298 ------------------ .../test_resources_manager.cpp | 29 +- 8 files changed, 298 insertions(+), 477 deletions(-) create mode 100644 src/amr/solvers/solver_field_evolvers.hpp delete mode 100644 src/amr/solvers/solver_ppc_model_view.hpp diff --git a/src/amr/multiphysics_integrator.hpp b/src/amr/multiphysics_integrator.hpp index 9c99bd3d5..80b8d163a 100644 --- a/src/amr/multiphysics_integrator.hpp +++ b/src/amr/multiphysics_integrator.hpp @@ -368,11 +368,6 @@ namespace solver } else load_balancer_manager_->estimate(*level, model); - - if (static_cast(levelNumber) == model_views_.size()) - model_views_.push_back(solver.make_view(*level, model)); - else - model_views_[levelNumber] = solver.make_view(*level, model); } @@ -391,9 +386,6 @@ namespace solver { messenger.registerLevel(hierarchy, ilvl); - model_views_.push_back(getSolver_(ilvl).make_view( - AMR_Types::getLevel(*hierarchy, ilvl), getModel_(ilvl))); - auto level = hierarchy->getPatchLevel(ilvl); for (auto& patch : *level) { @@ -533,8 +525,7 @@ namespace solver solver.prepareStep(model, *level, currentTime); fromCoarser.prepareStep(model, *level, currentTime); - solver.advanceLevel(*hierarchy, iLevel, getModelView_(iLevel), fromCoarser, currentTime, - newTime); + solver.advanceLevel(*hierarchy, iLevel, model, fromCoarser, currentTime, newTime); if (lastStep) { @@ -648,8 +639,6 @@ namespace solver std::vector>> solvers_; std::vector>> models_; - std::vector> model_views_; - std::vector> taggers_; std::map> messengers_; std::map> levelInitializers_; @@ -877,8 +866,6 @@ namespace solver } - auto& getModelView_(int iLevel) { return *model_views_[iLevel]; } - IPhysicalModel& getModel_(int iLevel) { diff --git a/src/amr/resources_manager/resources_manager.hpp b/src/amr/resources_manager/resources_manager.hpp index 12dfb8b2e..f660a2c92 100644 --- a/src/amr/resources_manager/resources_manager.hpp +++ b/src/amr/resources_manager/resources_manager.hpp @@ -343,6 +343,8 @@ namespace amr template struct LevelLooper { + using value_type = std::shared_ptr; + LevelLooper(ResourcesManager& rm, Level_t& lvl, Args&... arrgs) : rm{rm} , level{lvl} @@ -351,6 +353,23 @@ namespace amr } + struct Patch + { + Patch(LevelLooper* looper, std::shared_ptr const& patch) + : looper{looper} + , patch{patch} + { + looper->set(*patch); + } + + SAMRAI::hier::Patch const& operator*() { return *patch; } + + ~Patch() { looper->unset(); } + + LevelLooper* looper; + std::shared_ptr const patch; + }; + struct Iterator { void operator++() { ++raw; } @@ -380,6 +399,8 @@ namespace amr auto begin() { return Iterator{this, level.begin()}; } auto end() { return Iterator{this, level.end()}; }; + auto size() const { return static_cast(level.getLocalNumberOfPatches()); } + auto operator[](std::size_t const idx) { return Patch{this, level.getPatch(idx)}; } ResourcesManager& rm; Level_t& level; @@ -411,7 +432,8 @@ namespace amr if constexpr (has_compiletime_subresourcesview_list::value) { - // unpack the tuple subResources and apply for each element registerResources() + // unpack the tuple subResources and apply for each element + // registerResources() std::apply( [this, &IDs](auto&... subResource) { (this->getIDs_(subResource, IDs), ...); diff --git a/src/amr/solvers/solver.hpp b/src/amr/solvers/solver.hpp index 69624e735..8b43cf27e 100644 --- a/src/amr/solvers/solver.hpp +++ b/src/amr/solvers/solver.hpp @@ -13,20 +13,6 @@ #include -namespace PHARE::solver -{ - - -class ISolverModelView -{ -public: - using This = ISolverModelView; - - virtual ~ISolverModelView() = default; -}; - - -} // namespace PHARE::solver @@ -45,9 +31,10 @@ namespace solver class ISolver { public: - using patch_t = typename AMR_Types::patch_t; - using level_t = typename AMR_Types::level_t; - using hierarchy_t = typename AMR_Types::hierarchy_t; + using patch_t = AMR_Types::patch_t; + using level_t = AMR_Types::level_t; + using hierarchy_t = AMR_Types::hierarchy_t; + using IPhysicalModel_t = IPhysicalModel; /** * @brief return the name of the ISolver @@ -69,7 +56,7 @@ namespace solver * IPhysicalModel * @param model */ - virtual void registerResources(IPhysicalModel& model) = 0; + virtual void registerResources(IPhysicalModel_t& model) = 0; @@ -87,40 +74,36 @@ namespace solver * called before the advanceLevel() method. * */ - virtual void prepareStep(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, - double const currentTime) + virtual void prepareStep(IPhysicalModel_t& model, level_t& level, double const currentTime) = 0; /** * @brief accumulateFluxSum accumulates the flux sum(s) on the given PatchLevel for * refluxing later. */ - virtual void accumulateFluxSum(IPhysicalModel& model, - SAMRAI::hier::PatchLevel& level, double const coef) + virtual void accumulateFluxSum(IPhysicalModel_t& model, level_t& level, double const coef) = 0; /** * @brief resetFluxSum resets the flux sum(s) on the given PatchLevel to zero. */ - virtual void resetFluxSum(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level) - = 0; + virtual void resetFluxSum(IPhysicalModel_t& model, level_t& level) = 0; /** * @brief implements the reflux operations needed for a given solver. */ - virtual void reflux(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, - amr::IMessenger>& messenger, - double const time) + virtual void reflux(IPhysicalModel_t& model, level_t& level, + amr::IMessenger& messenger, double const time) = 0; /** * @brief advanceLevel advances the given level from t to t+dt */ virtual void advanceLevel(hierarchy_t const& hierarchy, int const levelNumber, - ISolverModelView& view, - amr::IMessenger>& fromCoarser, + IPhysicalModel_t& view, + amr::IMessenger& fromCoarser, double const currentTime, double const newTime) = 0; @@ -131,7 +114,7 @@ namespace solver * @brief allocate is used to allocate ISolver variables previously registered to the * ResourcesManager of the given model, onto the given Patch, at the given time. */ - virtual void allocate(IPhysicalModel& model, patch_t& patch, + virtual void allocate(IPhysicalModel_t& model, patch_t& patch, double const allocateTime) const = 0; @@ -143,9 +126,6 @@ namespace solver virtual ~ISolver() = default; - virtual std::shared_ptr make_view(level_t&, IPhysicalModel&) - = 0; - protected: explicit ISolver(std::string name) : solverName{std::move(name)} diff --git a/src/amr/solvers/solver_field_evolvers.hpp b/src/amr/solvers/solver_field_evolvers.hpp new file mode 100644 index 000000000..691111273 --- /dev/null +++ b/src/amr/solvers/solver_field_evolvers.hpp @@ -0,0 +1,117 @@ +#ifndef PHARE_AMR_SOLVERS_SOLVER_FIELD_EVOLVERS_HPP +#define PHARE_AMR_SOLVERS_SOLVER_FIELD_EVOLVERS_HPP + +#include "core/numerics/ohm/ohm.hpp" +#include "core/numerics/ampere/ampere.hpp" +#include "core/numerics/faraday/faraday.hpp" + +#include "amr/resources_manager/amr_utils.hpp" + + +namespace PHARE::solver +{ + + +template +class FaradayTransformer +{ + using core_type = PHARE::core::Faraday; + using level_t = AMR_Types::level_t; + +public: + void operator()(GridLayout& layout, auto&&... args) + { + auto _ = core::SetLayout(&layout, faraday_); + faraday_(args...); + } + + void operator()(level_t& level, auto& model, auto& B, auto& E, auto& Bnew, auto& dt) + { + auto& rm = *model.resourcesManager; + auto& state = model.state; + for (auto& patch : rm.enumerate(level, state, B, E, Bnew)) + { + auto layout = amr::layoutFromPatch(*patch); + (*this)(layout, B, E, Bnew, dt); + } + } + + + core_type faraday_; +}; + +template +class AmpereTransformer +{ + using core_type = PHARE::core::Ampere; + using level_t = AMR_Types::level_t; + +public: + void operator()(GridLayout& layout, auto&&... args) + { + auto _ = core::SetLayout(&layout, ampere_); + ampere_(args...); + } + + void operator()(level_t& level, auto& model, auto& B, auto& J) + { + auto& rm = *model.resourcesManager; + auto& state = model.state; + for (auto& patch : rm.enumerate(level, state, B, J)) + { + auto layout = amr::layoutFromPatch(*patch); + (*this)(layout, B, J); + } + } + + core_type ampere_; +}; + + +template +class OhmTransformer +{ + using core_type = PHARE::core::Ohm; + using level_t = AMR_Types::level_t; + +public: + explicit OhmTransformer(initializer::PHAREDict const& dict) + : ohm_{dict} + { + } + + void operator()(GridLayout& layout, auto&&... args) + { + auto _ = core::SetLayout(&layout, ohm_); + ohm_(args...); + } + + void operator()(level_t& level, auto& model, auto& B, auto& J, auto& E) + { + auto& electrons = model.state.electrons; + auto& rm = *model.resourcesManager; + for (auto& patch : rm.enumerate(level, electrons, B, J, E)) + { + auto layout = amr::layoutFromPatch(*patch); + auto& n = electrons.density(); + auto& Ve = electrons.velocity(); + auto& Pe = electrons.pressure(); + (*this)(layout, n, Ve, Pe, B, J, E); + } + } + + void operator()(level_t& level, auto& model, auto& B, auto& E) + { + (*this)(level, model, B, model.state.J, E); + } + + core_type ohm_; +}; + + + +} // namespace PHARE::solver + + + +#endif /* PHARE_AMR_SOLVERS_SOLVER_FIELD_EVOLVERS_HPP */ diff --git a/src/amr/solvers/solver_mhd.hpp b/src/amr/solvers/solver_mhd.hpp index e393ab725..615570772 100644 --- a/src/amr/solvers/solver_mhd.hpp +++ b/src/amr/solvers/solver_mhd.hpp @@ -14,9 +14,10 @@ namespace solver class SolverMHD : public ISolver { public: - using patch_t = typename AMR_Types::patch_t; - using level_t = typename AMR_Types::level_t; - using hierarchy_t = typename AMR_Types::hierarchy_t; + using patch_t = AMR_Types::patch_t; + using level_t = AMR_Types::level_t; + using hierarchy_t = AMR_Types::hierarchy_t; + using IPhysicalModel_t = IPhysicalModel; SolverMHD() : ISolver{"MHDSolver"} @@ -34,47 +35,38 @@ namespace solver } - void registerResources(IPhysicalModel& /*model*/) override {} + void registerResources(IPhysicalModel_t& /*model*/) override {} // TODO make this a resourcesUser - void allocate(IPhysicalModel& /*model*/, patch_t& /*patch*/, + void allocate(IPhysicalModel_t& /*model*/, patch_t& /*patch*/, double const /*allocateTime*/) const override { } - void prepareStep(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, + void prepareStep(IPhysicalModel_t& model, SAMRAI::hier::PatchLevel& level, double const currentTime) override { } - void accumulateFluxSum(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, + void accumulateFluxSum(IPhysicalModel_t& model, SAMRAI::hier::PatchLevel& level, double const coef) override { } - void resetFluxSum(IPhysicalModel& model, - SAMRAI::hier::PatchLevel& level) override - { - } + void resetFluxSum(IPhysicalModel_t& model, SAMRAI::hier::PatchLevel& level) override {} - virtual void reflux(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level, - amr::IMessenger>& messenger, + virtual void reflux(IPhysicalModel_t& model, SAMRAI::hier::PatchLevel& level, + amr::IMessenger& messenger, double const time) override { } void advanceLevel(hierarchy_t const& /*hierarchy*/, int const /*levelNumber*/, - ISolverModelView& /*view*/, - amr::IMessenger>& /*fromCoarser*/, + IPhysicalModel_t& /*view*/, + amr::IMessenger& /*fromCoarser*/, double const /*currentTime*/, double const /*newTime*/) override { } - - std::shared_ptr make_view(level_t&, IPhysicalModel&) override - { - throw std::runtime_error("Not implemented in mhd solver"); - return nullptr; - } }; } // namespace solver } // namespace PHARE diff --git a/src/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index 2779052fe..4bb516f78 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -1,6 +1,7 @@ #ifndef PHARE_SOLVER_PPC_HPP #define PHARE_SOLVER_PPC_HPP +#include "core/data/electrons/electrons.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep #include "core/numerics/ohm/ohm.hpp" @@ -14,7 +15,7 @@ #include "amr/solvers/solver.hpp" #include "amr/messengers/hybrid_messenger.hpp" #include "amr/resources_manager/amr_utils.hpp" -#include "amr/solvers/solver_ppc_model_view.hpp" +#include "amr/solvers/solver_field_evolvers.hpp" #include "amr/physical_models/physical_model.hpp" #include "amr/messengers/hybrid_messenger_info.hpp" @@ -46,10 +47,9 @@ class SolverPPC : public ISolver using IMessenger = amr::IMessenger; using HybridMessenger = amr::HybridMessenger; - using ModelViews_t = HybridPPCModelView; - using Faraday_t = ModelViews_t::Faraday_t; - using Ampere_t = ModelViews_t::Ampere_t; - using Ohm_t = ModelViews_t::Ohm_t; + using Faraday_t = FaradayTransformer; + using Ampere_t = AmpereTransformer; + using Ohm_t = OhmTransformer; using IonUpdater_t = PHARE::core::IonUpdater; Electromag electromagPred_{"EMPred"}; @@ -106,7 +106,7 @@ class SolverPPC : public ISolver void reflux(IPhysicalModel_t& model, SAMRAI::hier::PatchLevel& level, IMessenger& messenger, double const time) override; - void advanceLevel(hierarchy_t const& hierarchy, int const levelNumber, ISolverModelView& views, + void advanceLevel(hierarchy_t const& hierarchy, int const levelNumber, IPhysicalModel_t& views, IMessenger& fromCoarserMessenger, double const currentTime, double const newTime) override; @@ -118,10 +118,6 @@ class SolverPPC : public ISolver } - std::shared_ptr make_view(level_t& level, IPhysicalModel_t& model) override - { - return std::make_shared(level, dynamic_cast(model)); - } NO_DISCARD auto getCompileTimeResourcesViewList() { @@ -138,36 +134,36 @@ class SolverPPC : public ISolver using Messenger = amr::HybridMessenger; - void predictor1_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, + void predictor1_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime); - void predictor2_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, + void predictor2_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime); - void corrector_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, + void corrector_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime); - void average_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, - double const newTime); + void average_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const newTime); - void moveIons_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, + void moveIons_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime, core::UpdaterMode mode); struct TimeSetter { - template - void operator()(QuantityAccessor accessor) + void operator()(auto&... quantities) { - for (auto& state : views) - views.model().resourcesManager->setTime(accessor(state), *state.patch, newTime); + auto& rm = *model.resourcesManager; + for (auto& patch : rm.enumerate(level, quantities...)) + (model.resourcesManager->setTime(quantities, *patch, newTime), ...); } - ModelViews_t& views; + level_t& level; + HybridModel& model; double newTime; }; @@ -197,6 +193,13 @@ class SolverPPC : public ISolver return *level; } + void update_electrons(auto& level, auto& model) + { + auto& rm = *model.resourcesManager; + for (auto& patch : rm.enumerate(level, model.state.electrons)) + model.state.electrons.update(amr::layoutFromPatch(*patch)); + }; + using Boxing_t = core::UpdaterSelectionBoxing; std::unordered_map> boxing; @@ -335,16 +338,8 @@ void SolverPPC::reflux(IPhysicalModel_t& model, auto& hybridMessenger = dynamic_cast(messenger); auto& Eavg = electromagAvg_.E; auto& B = hybridModel.state.electromag.B; - - for (auto& patch : level) - { - core::Faraday faraday; - auto layout = amr::layoutFromPatch(*patch); - auto _sp = hybridModel.resourcesManager->setOnPatch(*patch, Bold_, Eavg, B); - auto _sl = core::SetLayout(&layout, faraday); - auto dt = time - oldTime_[level.getLevelNumber()]; - faraday(Bold_, Eavg, B, dt); - }; + auto dt = time - oldTime_[level.getLevelNumber()]; + faraday_(level, hybridModel, Bold_, Eavg, B, dt); hybridMessenger.fillMagneticGhosts(B, level, time); } @@ -353,102 +348,99 @@ void SolverPPC::reflux(IPhysicalModel_t& model, template void SolverPPC::advanceLevel(hierarchy_t const& hierarchy, - int const levelNumber, ISolverModelView& views, + int const levelNumber, + IPhysicalModel_t& level_model, IMessenger& fromCoarserMessenger, double const currentTime, double const newTime) { PHARE_LOG_SCOPE(1, "SolverPPC::advanceLevel"); - auto& modelView = dynamic_cast(views); + auto& model = dynamic_cast(level_model); auto& fromCoarser = dynamic_cast(fromCoarserMessenger); auto& level = setup_level(hierarchy, levelNumber); - predictor1_(level, modelView, fromCoarser, currentTime, newTime); + predictor1_(level, model, fromCoarser, currentTime, newTime); - average_(level, modelView, fromCoarser, newTime); + average_(level, model, fromCoarser, newTime); - moveIons_(level, modelView, fromCoarser, currentTime, newTime, core::UpdaterMode::domain_only); + moveIons_(level, model, fromCoarser, currentTime, newTime, core::UpdaterMode::domain_only); - predictor2_(level, modelView, fromCoarser, currentTime, newTime); + predictor2_(level, model, fromCoarser, currentTime, newTime); - average_(level, modelView, fromCoarser, newTime); + average_(level, model, fromCoarser, newTime); - moveIons_(level, modelView, fromCoarser, currentTime, newTime, core::UpdaterMode::all); + moveIons_(level, model, fromCoarser, currentTime, newTime, core::UpdaterMode::all); - corrector_(level, modelView, fromCoarser, currentTime, newTime); + corrector_(level, model, fromCoarser, currentTime, newTime); } template -void SolverPPC::predictor1_(level_t& level, ModelViews_t& views, +void SolverPPC::predictor1_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime) { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_"); - TimeSetter setTime{views, newTime}; + TimeSetter setTime{level, model, newTime}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.faraday"); auto dt = newTime - currentTime; - faraday_(views.layouts, views.electromag_B, views.electromag_E, views.electromagPred_B, dt); - setTime([](auto& state) -> auto& { return state.electromagPred.B; }); + faraday_(level, model, model.state.electromag.B, model.state.electromag.E, + electromagPred_.B, dt); + setTime(electromagPred_.B); fromCoarser.fillMagneticGhosts(electromagPred_.B, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.ampere"); - ampere_(views.layouts, views.electromagPred_B, views.J); - setTime([](auto& state) -> auto& { return state.J; }); - fromCoarser.fillCurrentGhosts(views.model().state.J, level, newTime); + ampere_(level, model, electromagPred_.B, model.state.J); + setTime(model.state.J); + fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.ohm"); - for (auto& state : views) - state.electrons.update(state.layout); - ohm_(views.layouts, views.N, views.Ve, views.Pe, views.electromagPred_B, views.J, - views.electromagPred_E); - setTime([](auto& state) -> auto& { return state.electromagPred.E; }); + update_electrons(level, model); + ohm_(level, model, electromagPred_.B, electromagPred_.E); + setTime(electromagPred_.E); } } template -void SolverPPC::predictor2_(level_t& level, ModelViews_t& views, +void SolverPPC::predictor2_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime) { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_"); - TimeSetter setTime{views, newTime}; + TimeSetter setTime{level, model, newTime}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.faraday"); auto dt = newTime - currentTime; - faraday_(views.layouts, views.electromag_B, views.electromagAvg_E, views.electromagPred_B, - dt); - setTime([](auto& state) -> auto& { return state.electromagPred.B; }); + faraday_(level, model, model.state.electromag.B, electromagAvg_.E, electromagPred_.B, dt); + setTime(electromagPred_.B); fromCoarser.fillMagneticGhosts(electromagPred_.B, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.ampere"); - ampere_(views.layouts, views.electromagPred_B, views.J); - setTime([](auto& state) -> auto& { return state.J; }); - fromCoarser.fillCurrentGhosts(views.model().state.J, level, newTime); + ampere_(level, model, electromagPred_.B, model.state.J); + setTime(model.state.J); + fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.ohm"); - for (auto& state : views) - state.electrons.update(state.layout); - ohm_(views.layouts, views.N, views.Ve, views.Pe, views.electromagPred_B, views.J, - views.electromagPred_E); - setTime([](auto& state) -> auto& { return state.electromagPred.E; }); + update_electrons(level, model); + ohm_(level, model, electromagPred_.B, electromagPred_.E); + setTime(electromagPred_.E); } } @@ -456,60 +448,61 @@ void SolverPPC::predictor2_(level_t& level, ModelViews_t template -void SolverPPC::corrector_(level_t& level, ModelViews_t& views, +void SolverPPC::corrector_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime) { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_"); auto levelNumber = level.getLevelNumber(); - TimeSetter setTime{views, newTime}; + TimeSetter setTime{level, model, newTime}; + auto& electromag = model.state.electromag; { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.faraday"); auto dt = newTime - currentTime; - faraday_(views.layouts, views.electromag_B, views.electromagAvg_E, views.electromag_B, dt); - setTime([](auto& state) -> auto& { return state.electromag.B; }); - fromCoarser.fillMagneticGhosts(views.model().state.electromag.B, level, newTime); + faraday_(level, model, electromag.B, electromagAvg_.E, electromag.B, dt); + setTime(model.state.electromag.B); + fromCoarser.fillMagneticGhosts(model.state.electromag.B, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.ampere"); - ampere_(views.layouts, views.electromag_B, views.J); - setTime([](auto& state) -> auto& { return state.J; }); - fromCoarser.fillCurrentGhosts(views.model().state.J, level, newTime); + ampere_(level, model, electromag.B, model.state.J); + setTime(model.state.J); + fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.ohm"); - for (auto& state : views) - state.electrons.update(state.layout); - ohm_(views.layouts, views.N, views.Ve, views.Pe, views.electromag_B, views.J, - views.electromag_E); - setTime([](auto& state) -> auto& { return state.electromag.E; }); - fromCoarser.fillElectricGhosts(views.model().state.electromag.E, level, newTime); + update_electrons(level, model); + ohm_(level, model, electromag.B, electromag.E); + setTime(model.state.electromag.E); + + fromCoarser.fillElectricGhosts(model.state.electromag.E, level, newTime); } } template -void SolverPPC::average_(level_t& level, ModelViews_t& views, +void SolverPPC::average_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const newTime) { PHARE_LOG_SCOPE(1, "SolverPPC::average_"); - TimeSetter setTime{views, newTime}; - - for (auto& state : views) + TimeSetter setTime{level, model, newTime}; + auto& rm = *model.resourcesManager; + auto& electromag = model.state.electromag; + for (auto& _ : rm.enumerate(level, electromag, electromagPred_, electromagAvg_)) { - PHARE::core::average(state.electromag.B, state.electromagPred.B, state.electromagAvg.B); - PHARE::core::average(state.electromag.E, state.electromagPred.E, state.electromagAvg.E); + PHARE::core::average(electromag.B, electromagPred_.B, electromagAvg_.B); + PHARE::core::average(electromag.E, electromagPred_.E, electromagAvg_.E); } - setTime([](auto& state) -> auto& { return state.electromagAvg.B; }); - setTime([](auto& state) -> auto& { return state.electromagAvg.E; }); + setTime(electromagAvg_.B); + setTime(electromagAvg_.E); // the following will fill E on all edges of all ghost cells, including those // on domain border. For level ghosts, electric field will be obtained from @@ -518,16 +511,15 @@ void SolverPPC::average_(level_t& level, ModelViews_t& v } -template -void _debug_log_move_ions(Args const&... args) +void _debug_log_move_ions(auto& level, auto& model) { - auto const& [views] = std::forward_as_tuple(args...); - + auto& rm = *model.resourcesManager; std::size_t nbrLevelGhostOldParticles = 0; std::size_t nbrLevelGhostParticles = 0; - for (auto& state : views) + + for (auto& patch : rm.enumerate(level, model.state.ions)) { - for (auto& pop : state.ions) + for (auto& pop : model.state.ions) { nbrLevelGhostOldParticles += pop.levelGhostParticlesOld().size(); nbrLevelGhostParticles += pop.levelGhostParticles().size(); @@ -544,23 +536,25 @@ void _debug_log_move_ions(Args const&... args) template -void SolverPPC::moveIons_(level_t& level, ModelViews_t& views, +void SolverPPC::moveIons_(level_t& level, HybridModel& model, Messenger& fromCoarser, double const currentTime, double const newTime, core::UpdaterMode mode) { PHARE_LOG_SCOPE(1, "SolverPPC::moveIons_"); - PHARE_DEBUG_DO(_debug_log_move_ions(views);) + PHARE_DEBUG_DO(_debug_log_move_ions(level, model);) - TimeSetter setTime{views, newTime}; + TimeSetter setTime{level, model, newTime}; + auto& rm = *model.resourcesManager; auto const& levelBoxing = boxing[level.getLevelNumber()]; + auto& ions = model.state.ions; try { auto dt = newTime - currentTime; - for (auto& state : views) - ionUpdater_.updatePopulations( - state.ions, state.electromagAvg, - levelBoxing.at(amr::to_string(state.patch->getGlobalId())), dt, mode); + for (auto& patch : rm.enumerate(level, ions, electromagAvg_)) + ionUpdater_.updatePopulations(ions, electromagAvg_, + levelBoxing.at(amr::to_string(patch->getGlobalId())), dt, + mode); } catch (core::DictionaryException const& ex) { @@ -570,22 +564,22 @@ void SolverPPC::moveIons_(level_t& level, ModelViews_t& throw core::DictionaryException{}("ID", "Updater::updatePopulations"); // this needs to be done before calling the messenger - setTime([](auto& state) -> auto& { return state.ions; }); + setTime(ions); - fromCoarser.fillFluxBorders(views.model().state.ions, level, newTime); - fromCoarser.fillDensityBorders(views.model().state.ions, level, newTime); - fromCoarser.fillIonPopMomentGhosts(views.model().state.ions, level, newTime); - fromCoarser.fillIonGhostParticles(views.model().state.ions, level, newTime); + fromCoarser.fillFluxBorders(model.state.ions, level, newTime); + fromCoarser.fillDensityBorders(model.state.ions, level, newTime); + fromCoarser.fillIonPopMomentGhosts(model.state.ions, level, newTime); + fromCoarser.fillIonGhostParticles(model.state.ions, level, newTime); - for (auto& state : views) - ionUpdater_.updateIons(state.ions); + for (auto& patch : rm.enumerate(level, ions)) + ionUpdater_.updateIons(ions); - fromCoarser.fillIonBorders(views.model().state.ions, level, newTime); + fromCoarser.fillIonBorders(model.state.ions, level, newTime); // no need to update time, since it has been done before // now Ni and Vi are calculated we can fill pure ghost nodes // these were not completed by the deposition of patch and levelghost particles - // fromCoarser.fillIonMomentGhosts(views.model().state.ions, level, newTime); + // fromCoarser.fillIonMomentGhosts(model.state.ions, level, newTime); } diff --git a/src/amr/solvers/solver_ppc_model_view.hpp b/src/amr/solvers/solver_ppc_model_view.hpp deleted file mode 100644 index c6c25e2ae..000000000 --- a/src/amr/solvers/solver_ppc_model_view.hpp +++ /dev/null @@ -1,298 +0,0 @@ -#ifndef PHARE_SOLVER_SOLVER_PPC_MODEL_VIEW_HPP -#define PHARE_SOLVER_SOLVER_PPC_MODEL_VIEW_HPP - -#include "core/numerics/ampere/ampere.hpp" -#include "core/numerics/faraday/faraday.hpp" -#include "core/numerics/ohm/ohm.hpp" - -#include "amr/solvers/solver.hpp" - -namespace PHARE::solver -{ -template -void assert_equal_sizes([[maybe_unused]] Vectors const&... vectors) -{ - PHARE_DEBUG_DO( // - auto tup = std::forward_as_tuple(vectors...); - assert(core::for_N_all>( - [&](auto i) { return std::get(tup).size() == std::get<0>(tup).size(); })); // - ) -} - -/*Faraday, Ampere, Ohm Transformers are abstraction that, from the solver viewpoint, act as Faraday, - * Ampere and Ohm algorithms, but take all patch views and hide the way these are processed, for - * instance to implement a parallelization decomposition*/ -template -class FaradayTransformer -{ - using core_type = PHARE::core::Faraday; - -public: - template - void operator()(GridLayouts const& layouts, VecFields const& B, VecFields const& E, - VecFields& Bnew, double dt) - { - assert_equal_sizes(B, E, Bnew); - for (std::size_t i = 0; i < B.size(); ++i) - { - auto _ = core::SetLayout(layouts[i], faraday_); - faraday_(*B[i], *E[i], *Bnew[i], dt); - } - } - - core_type faraday_; -}; - -template -class AmpereTransformer -{ - using core_type = PHARE::core::Ampere; - -public: - template - void operator()(GridLayouts const& layouts, VecFields const& B, VecFields& J) - { - assert_equal_sizes(B, J); - for (std::size_t i = 0; i < B.size(); ++i) - { - auto _ = core::SetLayout(layouts[i], ampere_); - ampere_(*B[i], *J[i]); - } - } - core_type ampere_; -}; - -template -class OhmTransformer -{ - using core_type = PHARE::core::Ohm; - -public: - explicit OhmTransformer(initializer::PHAREDict const& dict) - : ohm_{dict} - { - } - - template - void operator()(GridLayouts const& layouts, Fields const& n, VecFields const& Ve, - Fields const& Pe, VecFields const& B, VecFields const& J, VecFields& Enew) - { - assert_equal_sizes(n, Ve, Pe, B, J, Enew); - for (std::size_t i = 0; i < B.size(); ++i) - { - auto _ = core::SetLayout(layouts[i], ohm_); - ohm_(*n[i], *Ve[i], *Pe[i], *B[i], *J[i], *Enew[i]); - } - } - - core_type ohm_; -}; - - - -template -class HybridPPCModelView : public ISolverModelView -{ -public: - using This = HybridPPCModelView; - using HybridModel_t = HybridModel_; - using IPhysicalModel_t = HybridModel_t::Interface; - using patch_t = HybridModel_t::patch_t; - using Electrons = HybridModel_t::electrons_t; - using level_t = HybridModel_t::amr_types::level_t; - using Electromag = HybridModel_t::electromag_type; - using Ions = HybridModel_t::ions_type; - using ParticleArray_t = Ions::particle_array_type; - using Particle_t = ParticleArray_t::value_type; - using VecFieldT = HybridModel_t::vecfield_type; - using FieldT = HybridModel_t::field_type; - using GridLayout = HybridModel_t::gridlayout_type; - using Faraday_t = FaradayTransformer; - using Ampere_t = AmpereTransformer; - using Ohm_t = OhmTransformer; - - struct PatchState_t; - - using value_type = PatchState_t; - - template - struct iterator; - - HybridPPCModelView(level_t& level, IPhysicalModel_t& model) - : model_{dynamic_cast(model)} - { - onRegrid(level, model_); - } - - void onRegrid(level_t& level, HybridModel_t& hybridModel); - - auto begin() { return iterator{*this}; } - auto begin() const { return iterator{*this}; } - - auto end() { return iterator{*this, states.size()}; } - auto end() const { return iterator{*this, states.size()}; } - - auto& model() { return model_; } - auto& model() const { return model_; } - - HybridModel_t& model_; - std::vector> states; - - /*these vectors allow to access patch data per quantity - they are kept in synced with "states" so that client code (e.g. Faraday etc.) - can still have signatures revealing which quantities they use, and not take the states vector - entirely */ - std::vector electromag_E; - std::vector electromag_B; - std::vector electromagPred_E; - std::vector electromagPred_B; - std::vector electromagAvg_E; - std::vector electromagAvg_B; - std::vector J; - std::vector ions; - std::vector layouts; - std::vector Pe; - std::vector N; - std::vector Ve; - -private: - Electromag electromagPred_{"EMPred"}; - Electromag electromagAvg_{"EMAvg"}; - - auto vecs_tuple() - { - return std::forward_as_tuple(electromag_E, electromag_B, electromagPred_E, electromagPred_B, - electromagAvg_E, electromagAvg_B, Pe, N, Ve, J, ions, layouts); - } -}; - - -template -void HybridPPCModelView::onRegrid(level_t& level, HybridModel_t& hybridModel) -{ - auto& hybridState = hybridModel.state; - auto& rm = *hybridModel.resourcesManager; - - states.clear(); - - for (auto& patch : level) - { - { - auto _ = rm.setOnPatch(*patch, hybridState, electromagPred_, electromagAvg_); - states.emplace_back( // - PHARE::amr::layoutFromPatch(*patch), // - hybridState.ions, // - hybridState.J, // - hybridState.electromag, // - electromagPred_, // - electromagAvg_, // - hybridState.electrons, // - patch // - ); - } - assert(states.back().isUsable()); - } - core::apply(vecs_tuple(), [&](auto& vec) { vec.reserve(states.size()); }); - for (auto& state : states) - { - layouts.emplace_back(&state.layout); - ions.emplace_back(&state.ions); - J.emplace_back(&state.J); - Pe.emplace_back(&state.electrons.pressure()); - Ve.emplace_back(&state.electrons.velocity()); - N.emplace_back(&state.electrons.density()); - electromag_E.emplace_back(&state.electromag.E); - electromag_B.emplace_back(&state.electromag.B); - electromagAvg_E.emplace_back(&state.electromagAvg.E); - electromagAvg_B.emplace_back(&state.electromagAvg.B); - electromagPred_E.emplace_back(&state.electromagPred.E); - electromagPred_B.emplace_back(&state.electromagPred.B); - } -} - - -template -template -struct HybridPPCModelView::iterator -{ - using View = std::conditional_t, - HybridPPCModelView const>; - - iterator(View& view_, std::size_t const& i = 0) - : view{view_} - , idx{i} - { - } - - auto& operator*() - { - assert(view.states[idx].isUsable()); - return view.states[idx]; - } - auto& operator*() const - { - assert(view.states[idx].isUsable()); - return view.states[idx]; - } - - auto& operator++() - { - idx += 1; - return *this; - } - - bool operator!=(iterator const& that) const { return &view != &that.view or idx != that.idx; } - - View& view; - std::size_t idx = 0; -}; - - -template -struct HybridPPCModelView::PatchState_t -{ - GridLayout layout; - Ions ions; - VecFieldT J; - Electromag electromag, electromagPred, electromagAvg; - Electrons electrons; - std::shared_ptr patch; - - - bool isUsable() const - { - using Tuple = decltype((*this)()); - auto constexpr tuple_size = std::tuple_size_v; - Tuple tup = (*this)(); - - return core::for_N_all([&](auto i) { - auto const& [k, v] = std::get(tup); - auto isUsable = v->isUsable(); - if (!isUsable) - { - PHARE_LOG_LINE_STR("NOT USABLE: " << k); - } - return isUsable; - }); - } - - auto operator()() const - { - return core::make_named_tuple( // - std::make_pair("J", &J), // - std::make_pair("ions", &ions), // - std::make_pair("electrons", &electrons), // - std::make_pair("electromag", &electromag), // - std::make_pair("electromagAvg", &electromagAvg), // - std::make_pair("electromagPred", &electromagPred) // - ); - } -}; - - - -} // namespace PHARE::solver - - - -#endif /* PHARE_SOLVER_SOLVER_PPC_MODEL_VIEW_HPP */ diff --git a/tests/amr/resources_manager/test_resources_manager.cpp b/tests/amr/resources_manager/test_resources_manager.cpp index c202f9c99..cdac548ce 100644 --- a/tests/amr/resources_manager/test_resources_manager.cpp +++ b/tests/amr/resources_manager/test_resources_manager.cpp @@ -234,6 +234,33 @@ TYPED_TEST_P(aResourceUserCollection, hasPointersValidWithEnumerate) std::apply(check, resourceUserCollection); } +TYPED_TEST_P(aResourceUserCollection, hasPointersValidWithBracketOperator) +{ + TypeParam resourceUserCollection; + + auto check = [this](auto& resourceUserPack) { + auto& hierarchy_ = this->hierarchy->hierarchy; + auto& resourceUser = resourceUserPack.user; + auto& rm = this->resourcesManager; + for (int iLevel = 0; iLevel < hierarchy_->getNumberOfLevels(); ++iLevel) + { + auto patchLevel = hierarchy_->getPatchLevel(iLevel); + auto level_looper = rm.enumerate(*patchLevel, resourceUser); + + for (std::size_t i = 0; i < level_looper.size(); ++i) + { + auto const& _ = level_looper[i]; + EXPECT_TRUE(resourceUser.isUsable()); + EXPECT_FALSE(resourceUser.isSettable()); + } + EXPECT_FALSE(resourceUser.isUsable()); + EXPECT_TRUE(resourceUser.isSettable()); + } + }; + + std::apply(check, resourceUserCollection); +} + TEST(usingResourcesManager, toGetTimeOfAResourcesUser) { @@ -268,7 +295,7 @@ TEST(usingResourcesManager, toGetTimeOfAResourcesUser) REGISTER_TYPED_TEST_SUITE_P(aResourceUserCollection, hasPointersValidOnlyWithGuard, - hasPointersValidWithEnumerate); + hasPointersValidWithEnumerate, hasPointersValidWithBracketOperator); typedef ::testing::Types Date: Fri, 3 Apr 2026 13:39:56 +0200 Subject: [PATCH 2/3] drop LayoutHolder --- .../hybrid_level_initializer.hpp | 42 ++++------ src/amr/physical_models/hybrid_model.hpp | 15 ++-- src/amr/physical_models/physical_model.hpp | 28 +++++++ src/amr/solvers/solver_field_evolvers.hpp | 34 +++----- src/amr/solvers/solver_ppc.hpp | 16 ---- src/core/CMakeLists.txt | 1 - src/core/data/electrons/electrons.hpp | 9 +-- src/core/data/grid/gridlayout_utils.hpp | 48 ----------- src/core/numerics/ampere/ampere.hpp | 50 ++++++------ src/core/numerics/faraday/faraday.hpp | 43 +++++----- src/core/numerics/ohm/ohm.hpp | 61 +++++++------- tests/core/data/field/test_field_fixtures.hpp | 1 + .../data/vecfield/test_vecfield_fixtures.hpp | 2 + tests/core/numerics/ampere/test_main.cpp | 67 +++++----------- tests/core/numerics/faraday/test_main.cpp | 79 ++++--------------- tests/core/numerics/ohm/test_main.cpp | 65 ++++++--------- 16 files changed, 207 insertions(+), 354 deletions(-) delete mode 100644 src/core/data/grid/gridlayout_utils.hpp diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index 0249e3db6..2eb22896d 100644 --- a/src/amr/level_initializer/hybrid_level_initializer.hpp +++ b/src/amr/level_initializer/hybrid_level_initializer.hpp @@ -2,16 +2,14 @@ #define PHARE_HYBRID_LEVEL_INITIALIZER_HPP #include "core/errors.hpp" -#include "core/numerics/ohm/ohm.hpp" #include "core/utilities/mpi_utils.hpp" -#include "core/numerics/ampere/ampere.hpp" #include "core/numerics/moments/moments.hpp" -#include "core/data/grid/gridlayout_utils.hpp" #include "core/numerics/interpolator/interpolator.hpp" #include "amr/messengers/messenger.hpp" #include "amr/messengers/hybrid_messenger.hpp" #include "amr/resources_manager/amr_utils.hpp" +#include "amr/solvers/solver_field_evolvers.hpp" #include "amr/physical_models/physical_model.hpp" #include "amr/level_initializer/level_initializer.hpp" @@ -37,8 +35,11 @@ namespace solver static constexpr auto dimension = GridLayoutT::dimension; static constexpr auto interp_order = GridLayoutT::interp_order; - PHARE::core::Ohm ohm_; - PHARE::core::Ampere ampere_; + using Ampere_t = AmpereTransformer; + using Ohm_t = OhmTransformer; + + Ohm_t ohm_; + Ampere_t ampere_; inline bool isRootLevel(int const levelNumber) const { return levelNumber == 0; } @@ -47,6 +48,7 @@ namespace solver : ohm_{dict["algo"]["ohm"]} { } + virtual void initialize(std::shared_ptr const& hierarchy, int levelNumber, std::shared_ptr const& oldLevel, IPhysicalModelT& model, amr::IMessenger& messenger, double initDataTime, @@ -150,34 +152,22 @@ namespace solver if (!isRegriddingL0) if (isRootLevel(levelNumber)) { + TimeSetter setTime{level, hybridModel, 0.}; + auto& B = hybridModel.state.electromag.B; + auto& E = hybridModel.state.electromag.E; auto& J = hybridModel.state.J; - for (auto& patch : rm.enumerate(level, B, J)) - { - auto layout = PHARE::amr::layoutFromPatch(*patch); - auto __ = core::SetLayout(&layout, ampere_); - ampere_(B, J); - - rm.setTime(J, *patch, 0.); - } + ampere_(level, hybridModel, B, J); + setTime(J); hybMessenger.fillCurrentGhosts(J, level, 0.); auto& electrons = hybridModel.state.electrons; - auto& E = hybridModel.state.electromag.E; - - for (auto& patch : rm.enumerate(level, B, E, J, electrons)) - { - auto layout = PHARE::amr::layoutFromPatch(*patch); - electrons.update(layout); - auto& Ve = electrons.velocity(); - auto& Ne = electrons.density(); - auto& Pe = electrons.pressure(); - auto __ = core::SetLayout(&layout, ohm_); - ohm_(Ne, Ve, Pe, B, J, E); - rm.setTime(E, *patch, 0.); - } + for (auto& patch : rm.enumerate(level, electrons)) + electrons.update(amr::layoutFromPatch(*patch)); + ohm_(level, hybridModel, B, J, E); + setTime(E); hybMessenger.fillElectricGhosts(E, level, 0.); } diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index a485d0a1f..6f965ef4e 100644 --- a/src/amr/physical_models/hybrid_model.hpp +++ b/src/amr/physical_models/hybrid_model.hpp @@ -1,16 +1,19 @@ #ifndef PHARE_HYBRID_MODEL_HPP #define PHARE_HYBRID_MODEL_HPP -#include -#include "initializer/data_provider.hpp" +#include "core/def.hpp" #include "core/models/hybrid_state.hpp" -#include "amr/physical_models/physical_model.hpp" #include "core/data/ions/particle_initializers/particle_initializer_factory.hpp" -#include "amr/resources_manager/resources_manager.hpp" + +#include "initializer/data_provider.hpp" + +#include "amr/physical_models/physical_model.hpp" #include "amr/messengers/hybrid_messenger_info.hpp" -#include "core/data/vecfield/vecfield.hpp" -#include "core/def.hpp" +#include "amr/resources_manager/resources_manager.hpp" + + +#include namespace PHARE::solver { diff --git a/src/amr/physical_models/physical_model.hpp b/src/amr/physical_models/physical_model.hpp index cde754d02..f0bbdcfce 100644 --- a/src/amr/physical_models/physical_model.hpp +++ b/src/amr/physical_models/physical_model.hpp @@ -75,4 +75,32 @@ namespace solver } // namespace PHARE + + +namespace PHARE::solver +{ + + +template +struct TimeSetter +{ + void operator()(auto&... quantities) + { + auto& rm = *model.resourcesManager; + for (auto& patch : rm.enumerate(level, quantities...)) + (model.resourcesManager->setTime(quantities, *patch, newTime), ...); + } + + level_t& level; + Model& model; + double newTime; +}; + +template +TimeSetter(level_t&, Model&, double) -> TimeSetter; + + +} // namespace PHARE::solver + + #endif diff --git a/src/amr/solvers/solver_field_evolvers.hpp b/src/amr/solvers/solver_field_evolvers.hpp index 691111273..380281588 100644 --- a/src/amr/solvers/solver_field_evolvers.hpp +++ b/src/amr/solvers/solver_field_evolvers.hpp @@ -15,15 +15,11 @@ namespace PHARE::solver template class FaradayTransformer { - using core_type = PHARE::core::Faraday; + using core_type = core::Faraday; using level_t = AMR_Types::level_t; public: - void operator()(GridLayout& layout, auto&&... args) - { - auto _ = core::SetLayout(&layout, faraday_); - faraday_(args...); - } + void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } void operator()(level_t& level, auto& model, auto& B, auto& E, auto& Bnew, auto& dt) { @@ -35,23 +31,16 @@ class FaradayTransformer (*this)(layout, B, E, Bnew, dt); } } - - - core_type faraday_; }; template class AmpereTransformer { - using core_type = PHARE::core::Ampere; + using core_type = core::Ampere; using level_t = AMR_Types::level_t; public: - void operator()(GridLayout& layout, auto&&... args) - { - auto _ = core::SetLayout(&layout, ampere_); - ampere_(args...); - } + void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } void operator()(level_t& level, auto& model, auto& B, auto& J) { @@ -63,28 +52,23 @@ class AmpereTransformer (*this)(layout, B, J); } } - - core_type ampere_; }; template class OhmTransformer { - using core_type = PHARE::core::Ohm; + using info_type = core::OhmInfo; + using core_type = core::Ohm; using level_t = AMR_Types::level_t; public: explicit OhmTransformer(initializer::PHAREDict const& dict) - : ohm_{dict} + : info_{info_type::FROM(dict)} { } - void operator()(GridLayout& layout, auto&&... args) - { - auto _ = core::SetLayout(&layout, ohm_); - ohm_(args...); - } + void operator()(GridLayout& layout, auto&&... args) { core_type{info_, layout}(args...); } void operator()(level_t& level, auto& model, auto& B, auto& J, auto& E) { @@ -105,7 +89,7 @@ class OhmTransformer (*this)(level, model, B, model.state.J, E); } - core_type ohm_; + info_type info_; }; diff --git a/src/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index 4bb516f78..0b6b5cc25 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -9,7 +9,6 @@ #include "core/numerics/ampere/ampere.hpp" #include "core/data/vecfield/vecfield.hpp" #include "core/numerics/faraday/faraday.hpp" -#include "core/data/grid/gridlayout_utils.hpp" #include "core/numerics/ion_updater/ion_updater.hpp" #include "amr/solvers/solver.hpp" @@ -153,21 +152,6 @@ class SolverPPC : public ISolver double const currentTime, double const newTime, core::UpdaterMode mode); - struct TimeSetter - { - void operator()(auto&... quantities) - { - auto& rm = *model.resourcesManager; - for (auto& patch : rm.enumerate(level, quantities...)) - (model.resourcesManager->setTime(quantities, *patch, newTime), ...); - } - - level_t& level; - HybridModel& model; - double newTime; - }; - - void make_boxes(hierarchy_t const& hierarchy, level_t& level) { int const lvlNbr = level.getLevelNumber(); diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 10273d690..5494f416f 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -8,7 +8,6 @@ set( SOURCES_INC data/grid/gridlayout.hpp data/grid/gridlayout_impl.hpp data/grid/gridlayoutimplyee.hpp - data/grid/gridlayout_utils.hpp data/ndarray/ndarray_vector.hpp data/particles/particle.hpp data/particles/particle_utilities.hpp diff --git a/src/core/data/electrons/electrons.hpp b/src/core/data/electrons/electrons.hpp index b18d0c290..62d817813 100644 --- a/src/core/data/electrons/electrons.hpp +++ b/src/core/data/electrons/electrons.hpp @@ -1,16 +1,11 @@ #ifndef PHARE_ELECTRONS_HPP #define PHARE_ELECTRONS_HPP +#include "core/def.hpp" #include "core/hybrid/hybrid_quantities.hpp" #include "core/data/vecfield/vecfield_component.hpp" -#include "core/data/grid/gridlayout_utils.hpp" -#include "core/data/grid/gridlayoutdefs.hpp" -#include "core/utilities/index/index.hpp" -#include "core/def.hpp" -#include "core/logger.hpp" #include "initializer/data_provider.hpp" -#include namespace PHARE::core @@ -293,7 +288,7 @@ class ElectronMomentModel template -class Electrons : public LayoutHolder +class Electrons { using VecField = typename Ions::vecfield_type; using Field = typename Ions::field_type; diff --git a/src/core/data/grid/gridlayout_utils.hpp b/src/core/data/grid/gridlayout_utils.hpp deleted file mode 100644 index 3574acb99..000000000 --- a/src/core/data/grid/gridlayout_utils.hpp +++ /dev/null @@ -1,48 +0,0 @@ -#ifndef PHARE_GRIDLAYOUT_UTILS_HPP -#define PHARE_GRIDLAYOUT_UTILS_HPP - -#include "core/def.hpp" -#include -#include - -namespace PHARE::core -{ -template -class LayoutHolder -{ -protected: - GridLayout* layout_{nullptr}; - -public: - void setLayout(GridLayout* ptr) { layout_ = ptr; } - - NO_DISCARD bool hasLayout() const { return layout_ != nullptr; } -}; - - - - -template -class SetLayout -{ -public: - SetLayout(GridLayout* ptr, GridLayoutSettable&... settables) - : settables_{settables...} - { - std::apply([ptr](auto&... settable) { (settable.setLayout(ptr), ...); }, settables_); - } - - ~SetLayout() - { - std::apply([](auto&... settable) { (settable.setLayout(nullptr), ...); }, settables_); - } - - -private: - std::tuple settables_; -}; -} // namespace PHARE::core - - - -#endif diff --git a/src/core/numerics/ampere/ampere.hpp b/src/core/numerics/ampere/ampere.hpp index 1f973d011..d3c986a20 100644 --- a/src/core/numerics/ampere/ampere.hpp +++ b/src/core/numerics/ampere/ampere.hpp @@ -1,44 +1,46 @@ #ifndef PHARE_CORE_NUMERICS_AMPERE_AMPERE_HPP #define PHARE_CORE_NUMERICS_AMPERE_AMPERE_HPP -#include -#include - #include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/grid/gridlayout_utils.hpp" #include "core/data/vecfield/vecfield_component.hpp" -#include "core/utilities/index/index.hpp" - namespace PHARE::core { + + + + template -class Ampere : public LayoutHolder +class Ampere { constexpr static auto dimension = GridLayout::dimension; - using LayoutHolder::layout_; public: + Ampere(GridLayout const& layout) + : layout_{layout} + { + } + + template void operator()(VecField const& B, VecField& J) { - if (!this->hasLayout()) - throw std::runtime_error( - "Error - Ampere - GridLayout not set, cannot proceed to calculate ampere()"); - // can't use structured bindings because // "reference to local binding declared in enclosing function" auto& Jx = J(Component::X); auto& Jy = J(Component::Y); auto& Jz = J(Component::Z); - layout_->evalOnBox(Jx, [&](auto&... args) mutable { JxEq_(Jx, B, args...); }); - layout_->evalOnBox(Jy, [&](auto&... args) mutable { JyEq_(Jy, B, args...); }); - layout_->evalOnBox(Jz, [&](auto&... args) mutable { JzEq_(Jz, B, args...); }); + layout_.evalOnBox(Jx, [&](auto&... args) mutable { JxEq_(Jx, B, args...); }); + layout_.evalOnBox(Jy, [&](auto&... args) mutable { JyEq_(Jy, B, args...); }); + layout_.evalOnBox(Jz, [&](auto&... args) mutable { JzEq_(Jz, B, args...); }); } private: + GridLayout layout_; + + template void JxEq_(Field& Jx, VecField const& B, Indexes const&... ijk) const { @@ -48,11 +50,11 @@ class Ampere : public LayoutHolder Jx(ijk...) = 0.0; if constexpr (dimension == 2) - Jx(ijk...) = layout_->template deriv(Bz, {ijk...}); + Jx(ijk...) = layout_.template deriv(Bz, {ijk...}); if constexpr (dimension == 3) - Jx(ijk...) = layout_->template deriv(Bz, {ijk...}) - - layout_->template deriv(By, {ijk...}); + Jx(ijk...) = layout_.template deriv(Bz, {ijk...}) + - layout_.template deriv(By, {ijk...}); } template @@ -61,11 +63,11 @@ class Ampere : public LayoutHolder auto const& [Bx, By, Bz] = B(); if constexpr (dimension == 1 || dimension == 2) - Jy(ijk...) = -layout_->template deriv(Bz, {ijk...}); + Jy(ijk...) = -layout_.template deriv(Bz, {ijk...}); if constexpr (dimension == 3) - Jy(ijk...) = layout_->template deriv(Bx, {ijk...}) - - layout_->template deriv(Bz, {ijk...}); + Jy(ijk...) = layout_.template deriv(Bx, {ijk...}) + - layout_.template deriv(Bz, {ijk...}); } template @@ -74,11 +76,11 @@ class Ampere : public LayoutHolder auto const& [Bx, By, Bz] = B(); if constexpr (dimension == 1) - Jz(ijk...) = layout_->template deriv(By, {ijk...}); + Jz(ijk...) = layout_.template deriv(By, {ijk...}); else - Jz(ijk...) = layout_->template deriv(By, {ijk...}) - - layout_->template deriv(Bx, {ijk...}); + Jz(ijk...) = layout_.template deriv(By, {ijk...}) + - layout_.template deriv(Bx, {ijk...}); } }; diff --git a/src/core/numerics/faraday/faraday.hpp b/src/core/numerics/faraday/faraday.hpp index 0fe28cc77..557e3222d 100644 --- a/src/core/numerics/faraday/faraday.hpp +++ b/src/core/numerics/faraday/faraday.hpp @@ -4,26 +4,29 @@ #include #include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/grid/gridlayout_utils.hpp" #include "core/data/vecfield/vecfield_component.hpp" namespace PHARE::core { + + + template -class Faraday : public LayoutHolder +class Faraday { constexpr static auto dimension = GridLayout::dimension; - using LayoutHolder::layout_; public: + Faraday(GridLayout const& layout) + : layout_{layout} + { + } + + template void operator()(VecField const& B, VecField const& E, VecField& Bnew, double dt) { - if (!this->hasLayout()) - throw std::runtime_error( - "Error - Faraday - GridLayout not set, cannot proceed to calculate faraday()"); - if (!(B.isUsable() && E.isUsable() && Bnew.isUsable())) throw std::runtime_error("Error - Faraday - not all VecField parameters are usable"); @@ -39,15 +42,15 @@ class Faraday : public LayoutHolder auto& Bynew = Bnew(Component::Y); auto& Bznew = Bnew(Component::Z); - layout_->evalOnBox(Bxnew, [&](auto&... args) mutable { BxEq_(Bx, E, Bxnew, args...); }); - layout_->evalOnBox(Bynew, [&](auto&... args) mutable { ByEq_(By, E, Bynew, args...); }); - layout_->evalOnBox(Bznew, [&](auto&... args) mutable { BzEq_(Bz, E, Bznew, args...); }); + layout_.evalOnBox(Bxnew, [&](auto&... args) mutable { BxEq_(Bx, E, Bxnew, args...); }); + layout_.evalOnBox(Bynew, [&](auto&... args) mutable { ByEq_(By, E, Bynew, args...); }); + layout_.evalOnBox(Bznew, [&](auto&... args) mutable { BzEq_(Bz, E, Bznew, args...); }); } private: double dt_; - + GridLayout layout_; template void BxEq_(Field const& Bx, VecField const& E, Field& Bxnew, Indexes const&... ijk) const @@ -58,11 +61,11 @@ class Faraday : public LayoutHolder Bxnew(ijk...) = Bx(ijk...); if constexpr (dimension == 2) - Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_->template deriv(Ez, {ijk...}); + Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_.template deriv(Ez, {ijk...}); if constexpr (dimension == 3) - Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_->template deriv(Ez, {ijk...}) - + dt_ * layout_->template deriv(Ey, {ijk...}); + Bxnew(ijk...) = Bx(ijk...) - dt_ * layout_.template deriv(Ez, {ijk...}) + + dt_ * layout_.template deriv(Ey, {ijk...}); } template @@ -71,11 +74,11 @@ class Faraday : public LayoutHolder auto const& [Ex, _, Ez] = E(); if constexpr (dimension == 1 || dimension == 2) - Bynew(ijk...) = By(ijk...) + dt_ * layout_->template deriv(Ez, {ijk...}); + Bynew(ijk...) = By(ijk...) + dt_ * layout_.template deriv(Ez, {ijk...}); if constexpr (dimension == 3) - Bynew(ijk...) = By(ijk...) - dt_ * layout_->template deriv(Ex, {ijk...}) - + dt_ * layout_->template deriv(Ez, {ijk...}); + Bynew(ijk...) = By(ijk...) - dt_ * layout_.template deriv(Ex, {ijk...}) + + dt_ * layout_.template deriv(Ez, {ijk...}); } template @@ -84,11 +87,11 @@ class Faraday : public LayoutHolder auto const& [Ex, Ey, _] = E(); if constexpr (dimension == 1) - Bznew(ijk...) = Bz(ijk...) - dt_ * layout_->template deriv(Ey, {ijk...}); + Bznew(ijk...) = Bz(ijk...) - dt_ * layout_.template deriv(Ey, {ijk...}); else - Bznew(ijk...) = Bz(ijk...) - dt_ * layout_->template deriv(Ey, {ijk...}) - + dt_ * layout_->template deriv(Ex, {ijk...}); + Bznew(ijk...) = Bz(ijk...) - dt_ * layout_.template deriv(Ey, {ijk...}) + + dt_ * layout_.template deriv(Ex, {ijk...}); } }; diff --git a/src/core/numerics/ohm/ohm.hpp b/src/core/numerics/ohm/ohm.hpp index 107c91a8e..856558872 100644 --- a/src/core/numerics/ohm/ohm.hpp +++ b/src/core/numerics/ohm/ohm.hpp @@ -4,30 +4,42 @@ #include "core/utilities/index/index.hpp" #include "core/data/grid/gridlayoutdefs.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 }; +struct OhmInfo +{ + double const eta_; + double const nu_; + HyperMode const hyper_mode; + + OhmInfo static FROM(initializer::PHAREDict const& dict) + { + return {dict["resistivity"].template to(), + dict["hyper_resistivity"].template to(), + cppdict::get_value(dict, "hyper_mode", std::string{"constant"}) == "constant" + ? HyperMode::constant + : HyperMode::spatial}; + } +}; + template -class Ohm : public LayoutHolder +class Ohm : public OhmInfo { + using Super = OhmInfo; constexpr static auto dimension = GridLayout::dimension; - using LayoutHolder::layout_; public: - 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} + explicit Ohm(OhmInfo const& info, GridLayout const& layout) + : Super{info} + , layout_{layout} { } @@ -37,30 +49,21 @@ class Ohm : public LayoutHolder { using Pack = OhmPack; - if (!this->hasLayout()) - throw std::runtime_error( - "Error - Ohm - GridLayout not set, cannot proceed to calculate ohm()"); - auto const& [Exnew, Eynew, Eznew] = Enew(); - layout_->evalOnBox(Exnew, [&](auto&... args) mutable { + layout_.evalOnBox(Exnew, [&](auto&... args) mutable { this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); }); - layout_->evalOnBox(Eynew, [&](auto&... args) mutable { + layout_.evalOnBox(Eynew, [&](auto&... args) mutable { this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); }); - layout_->evalOnBox(Eznew, [&](auto&... args) mutable { + layout_.evalOnBox(Eznew, [&](auto&... args) mutable { this->template E_Eq_(Pack{Enew, n, Pe, Ve, B, J}, args...); }); } - - - private: - double const eta_; - double const nu_; - HyperMode const hyper_mode; + GridLayout layout_; template struct OhmPack @@ -270,7 +273,7 @@ class Ohm : public LayoutHolder { auto const nOnEx = GridLayout::project(n, index, GridLayout::momentsToEx()); - auto gradPOnEx = layout_->template deriv(Pe, index); // TODO : issue 3391 + auto gradPOnEx = layout_.template deriv(Pe, index); // TODO : issue 3391 return -gradPOnEx / nOnEx; } @@ -282,7 +285,7 @@ class Ohm : public LayoutHolder auto const nOnEy = GridLayout::project(n, index, GridLayout::momentsToEy()); auto gradPOnEy - = layout_->template deriv(Pe, index); // TODO : issue 3391 + = layout_.template deriv(Pe, index); // TODO : issue 3391 return -gradPOnEy / nOnEy; } @@ -299,7 +302,7 @@ class Ohm : public LayoutHolder auto const nOnEz = GridLayout::project(n, index, GridLayout::momentsToEz()); auto gradPOnEz - = layout_->template deriv(Pe, index); // TODO : issue 3391 + = layout_.template deriv(Pe, index); // TODO : issue 3391 return -gradPOnEz / nOnEz; } @@ -353,7 +356,7 @@ class Ohm : public LayoutHolder 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); + return -nu_ * layout_.laplacian(J(component), index); } @@ -361,7 +364,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 = 1. / std::pow(4, layout_->levelNumber()); + auto const lvlCoeff = 1. / std::pow(4, layout_.levelNumber()); auto constexpr min_density = 0.1; auto computeHR = [&](auto BxProj, auto ByProj, auto BzProj, auto nProj) { @@ -371,7 +374,7 @@ class Ohm : public LayoutHolder auto const nOnE = GridLayout::project(n, index, nProj); auto b = std::sqrt(BxOnE * BxOnE + ByOnE * ByOnE + BzOnE * BzOnE); return -nu_ * (b / (nOnE + min_density) + 1) * lvlCoeff - * layout_->laplacian(J(component), index); + * layout_.laplacian(J(component), index); }; if constexpr (component == Component::X) diff --git a/tests/core/data/field/test_field_fixtures.hpp b/tests/core/data/field/test_field_fixtures.hpp index c63c1ce97..60fabbbec 100644 --- a/tests/core/data/field/test_field_fixtures.hpp +++ b/tests/core/data/field/test_field_fixtures.hpp @@ -2,6 +2,7 @@ #define PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP #include "core/data/field/field.hpp" +#include "core/hybrid/hybrid_quantities.hpp" namespace PHARE::core { diff --git a/tests/core/data/vecfield/test_vecfield_fixtures.hpp b/tests/core/data/vecfield/test_vecfield_fixtures.hpp index 9caa093c7..257243ec1 100644 --- a/tests/core/data/vecfield/test_vecfield_fixtures.hpp +++ b/tests/core/data/vecfield/test_vecfield_fixtures.hpp @@ -1,6 +1,8 @@ #ifndef PHARE_TEST_CORE_DATA_TEST_VECFIELD_FIXTURES_HPP #define PHARE_TEST_CORE_DATA_TEST_VECFIELD_FIXTURES_HPP +#include "core/data/vecfield/vecfield.hpp" + #include "tests/core/data/field/test_field_fixtures.hpp" #include "tests/core/data/tensorfield/test_tensorfield_fixtures.hpp" diff --git a/tests/core/numerics/ampere/test_main.cpp b/tests/core/numerics/ampere/test_main.cpp index 738c8f588..ba5747a6a 100644 --- a/tests/core/numerics/ampere/test_main.cpp +++ b/tests/core/numerics/ampere/test_main.cpp @@ -1,24 +1,19 @@ -#include "gmock/gmock.h" -#include "gtest/gtest.h" - -#include -#include -#include "core/data/grid/grid.hpp" #include "core/data/field/field.hpp" #include "core/data/grid/gridlayout.hpp" -#include "core/data/grid/gridlayout_impl.hpp" -#include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/vecfield/vecfield.hpp" -#include "core/numerics/ampere/ampere.hpp" -#include "core/utilities/box/box.hpp" #include "core/utilities/index/index.hpp" +#include "core/numerics/ampere/ampere.hpp" +#include "core/data/grid/gridlayoutdefs.hpp" +#include "core/data/grid/gridlayoutimplyee.hpp" #include "tests/core/data/field/test_field.hpp" -#include "tests/core/data/vecfield/test_vecfield.hpp" #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" -#include "tests/core/data/gridlayout/gridlayout_test.hpp" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include using namespace PHARE::core; @@ -73,38 +68,19 @@ struct GridLayoutMock3D TEST(Ampere, canBe1D) { - Ampere ampere; + Ampere ampere{GridLayoutMock1D{}}; } TEST(Ampere, canBe2D) { - Ampere ampere; + Ampere ampere{GridLayoutMock2D{}}; } TEST(Ampere, canBe3D) { - Ampere ampere; + Ampere ampere{GridLayoutMock3D{}}; } -TEST(Ampere, shouldBeGivenAGridLayoutPointerToBeOperational) -{ - std::size_t constexpr interp = 1; - - auto constexpr check = [](auto ic) { - auto constexpr dim = ic(); - - VecFieldMock> B, J; - using GridLayout = GridLayout>; - Ampere ampere; - auto layout = std::make_unique>(); - EXPECT_ANY_THROW(ampere(B, J)); - ampere.setLayout(layout.get()); - }; - - check(std::integral_constant{}); - check(std::integral_constant{}); - check(std::integral_constant{}); -} std::vector read(std::string filename) @@ -127,12 +103,12 @@ class Ampere1DTest : public ::testing::Test static constexpr std::size_t interp_order = 1; using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; + using Ampere_t = Ampere>; + GridLayout layout; UsableVecFieldND B, J; - Ampere> ampere; - public: Ampere1DTest() : layout{{{0.1}}, {{50}}, Point{0.}} @@ -150,12 +126,12 @@ class Ampere2DTest : public ::testing::Test static constexpr std::size_t interp_order = 1; using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; + using Ampere_t = Ampere>; + GridLayout layout; UsableVecFieldND B, J; - Ampere> ampere; - public: Ampere2DTest() : layout{{{0.1, 0.2}}, {{50, 30}}, Point{0., 0.}} @@ -173,12 +149,12 @@ class Ampere3DTest : public ::testing::Test static constexpr std::size_t interp_order = 1; using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; + using Ampere_t = Ampere>; + GridLayout layout; UsableVecFieldND B, J; - Ampere> ampere; - public: Ampere3DTest() : layout{{{0.1, 0.2, 0.3}}, {{50, 30, 40}}, Point{0., 0., 0.}} @@ -211,8 +187,7 @@ TEST_F(Ampere1DTest, ampere1DCalculatedOk) Bz(ix) = std::sin(2 * M_PI / 5. * point[0]); } - ampere.setLayout(&layout); - ampere(B, J); + Ampere_t{layout}(B, J); auto psi_p_X = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto pei_p_X = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); @@ -276,8 +251,7 @@ TEST_F(Ampere2DTest, ampere2DCalculatedOk) } } - ampere.setLayout(&layout); - ampere(B, J); + Ampere_t{layout}(B, J); auto psi_p_X = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto pei_p_X = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); @@ -394,8 +368,7 @@ TEST_F(Ampere3DTest, ampere3DCalculatedOk) } } - ampere.setLayout(&layout); - ampere(B, J); + Ampere_t{layout}(B, J); auto psi_p_X = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto pei_p_X = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); diff --git a/tests/core/numerics/faraday/test_main.cpp b/tests/core/numerics/faraday/test_main.cpp index 8c437ba8a..30f3acf12 100644 --- a/tests/core/numerics/faraday/test_main.cpp +++ b/tests/core/numerics/faraday/test_main.cpp @@ -1,28 +1,22 @@ -#include "gmock/gmock.h" -#include "gtest/gtest.h" - -#include -#include -#include "core/data/grid/grid.hpp" #include "core/data/grid/gridlayout.hpp" -#include "core/data/grid/gridlayout_impl.hpp" -#include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/vecfield/vecfield.hpp" -#include "core/numerics/faraday/faraday.hpp" -#include "core/utilities/box/box.hpp" #include "core/utilities/index/index.hpp" #include "core/utilities/point/point.hpp" +#include "core/data/grid/gridlayoutdefs.hpp" +#include "core/numerics/faraday/faraday.hpp" +#include "core/data/grid/gridlayoutimplyee.hpp" #include "tests/core/data/field/test_field.hpp" -#include "tests/core/data/vecfield/test_vecfield.hpp" #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" -#include "tests/core/data/gridlayout/gridlayout_test.hpp" -using namespace PHARE::core; +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +using namespace PHARE::core; struct GridLayoutMock1D @@ -72,54 +66,22 @@ struct GridLayoutMock3D TEST(Faraday, canBe1D) { - Faraday faraday; + Faraday faraday{GridLayoutMock1D{}}; } TEST(Faraday, canBe2D) { - Faraday faraday; + Faraday faraday{GridLayoutMock2D{}}; } TEST(Faraday, canBe3D) { - Faraday faraday; + Faraday faraday{GridLayoutMock3D{}}; } -TEST(Faraday, shouldBeGivenAGridLayoutPointerToBeOperational) -{ - { - using GridLayout = GridLayout>; - VecFieldMock> B_1, E_1, Bnew_1; - Faraday faraday1d; - auto layout1d = std::make_unique>(); - EXPECT_ANY_THROW(faraday1d(B_1, E_1, Bnew_1, 1.)); - faraday1d.setLayout(layout1d.get()); - } - - { - using GridLayout = GridLayout>; - VecFieldMock> B_2, E_2, Bnew_2; - Faraday faraday2d; - auto layout2d = std::make_unique>(); - EXPECT_ANY_THROW(faraday2d(B_2, E_2, Bnew_2, 1.)); - faraday2d.setLayout(layout2d.get()); - } - - { - using GridLayout = GridLayout>; - VecFieldMock> B_3, E_3, Bnew_3; - Faraday faraday3d; - auto layout3d = std::make_unique>(); - EXPECT_ANY_THROW(faraday3d(B_3, E_3, Bnew_3, 1.)); - faraday3d.setLayout(layout3d.get()); - } -} - - - std::vector read(std::string filename) { @@ -142,13 +104,12 @@ class Faraday1DTest : public ::testing::Test using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; + using Faraday_t = Faraday>; GridLayout layout; UsableVecFieldND B, E, Bnew; - Faraday> faraday; - public: Faraday1DTest() : layout{{{0.1}}, {{50}}, Point{0.}} @@ -170,13 +131,12 @@ class Faraday2DTest : public ::testing::Test using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; + using Faraday_t = Faraday>; GridLayout layout; UsableVecFieldND B, E, Bnew; - Faraday> faraday; - public: Faraday2DTest() : layout{{{0.1, 0.2}}, {{50, 30}}, Point{0., 0.}} @@ -198,13 +158,11 @@ class Faraday3DTest : public ::testing::Test using UsableVecFieldND = UsableVecField; using GridLayoutImpl = GridLayoutImplYee; - + using Faraday_t = Faraday>; GridLayout layout; UsableVecFieldND B, E, Bnew; - Faraday> faraday; - public: Faraday3DTest() : layout{{{0.1, 0.2, 0.3}}, {{50, 30, 40}}, Point{0., 0., 0.}} @@ -253,8 +211,7 @@ TEST_F(Faraday1DTest, Faraday1DCalculatedOk) Bz(ix) = std::tanh(point[0] - 5. / 2.); } - faraday.setLayout(&layout); - faraday(B, E, Bnew, 1.); + Faraday_t{layout}(B, E, Bnew, 1.); auto psi_d_X = this->layout.physicalStartIndex(QtyCentering::dual, Direction::X); auto pei_d_X = this->layout.physicalEndIndex(QtyCentering::dual, Direction::X); @@ -358,8 +315,7 @@ TEST_F(Faraday2DTest, Faraday2DCalculatedOk) } } - faraday.setLayout(&layout); - faraday(B, E, Bnew, 1.); + Faraday_t{layout}(B, E, Bnew, 1.); auto psi_p_X = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto pei_p_X = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); @@ -528,8 +484,7 @@ TEST_F(Faraday3DTest, Faraday3DCalculatedOk) } } - faraday.setLayout(&layout); - faraday(B, E, Bnew, 1.); + Faraday_t{layout}(B, E, Bnew, 1.); auto psi_p_X = this->layout.physicalStartIndex(QtyCentering::primal, Direction::X); auto pei_p_X = this->layout.physicalEndIndex(QtyCentering::primal, Direction::X); diff --git a/tests/core/numerics/ohm/test_main.cpp b/tests/core/numerics/ohm/test_main.cpp index 7ca7c1718..4dbba10c5 100644 --- a/tests/core/numerics/ohm/test_main.cpp +++ b/tests/core/numerics/ohm/test_main.cpp @@ -1,22 +1,23 @@ -#include "gmock/gmock.h" -#include "gtest/gtest.h" - -#include -#include - -#include "core/data/field/field.hpp" +// #include "phare_core.hpp" +// #include "core/data/field/field.hpp" #include "core/data/grid/gridlayout.hpp" #include "core/data/grid/gridlayout_impl.hpp" #include "core/data/grid/gridlayoutdefs.hpp" -#include "core/data/vecfield/vecfield.hpp" +// #include "core/data/vecfield/vecfield.hpp" #include "core/numerics/ohm/ohm.hpp" -#include "core/utilities/index/index.hpp" +// #include "core/utilities/index/index.hpp" -#include "phare_core.hpp" #include "tests/core/data/vecfield/test_vecfield_fixtures.hpp" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#include +#include + using namespace PHARE::core; @@ -69,13 +70,15 @@ struct OhmTest : public ::testing::Test using GridYee = GridLayout>; using UsableVecFieldND = UsableVecField; using Grid_t = Grid, HybridQuantity::Scalar>; - GridYee layout = NDlayout::create(); + using Ohm_t = Ohm; + + GridYee layout = NDlayout::create(); Grid_t n; Grid_t P; UsableVecFieldND V, B, J, Enew; - Ohm ohm; + OhmInfo ohm_info; OhmTest() : n{"n", HybridQuantity::Scalar::rho, layout.allocSize(HybridQuantity::Scalar::rho)} @@ -84,7 +87,7 @@ struct OhmTest : public ::testing::Test , B{"B", layout, HybridQuantity::Vector::B} , J{"J", layout, HybridQuantity::Vector::J} , Enew{"Enew", layout, HybridQuantity::Vector::E} - , ohm{createDict()} + , ohm_info{OhmInfo::FROM(createDict())} { auto const& [Bx, By, Bz] = B(); auto const& [Jx, Jy, Jz] = J(); @@ -292,37 +295,13 @@ TYPED_TEST_SUITE(OhmTest, OhmTupleInfos); -TYPED_TEST(OhmTest, ThatOhmHasCtorWithDict) +TYPED_TEST(OhmTest, ThatOhmInfoCanBeBuiltFromDict) { - TypeParam pair; - auto constexpr dim = pair.first(); - auto constexpr interp = pair.second(); - - using GridYee = GridLayout>; - - Ohm ohm(createDict()); + auto info = OhmInfo::FROM(createDict()); } - -TYPED_TEST(OhmTest, ShouldBeGivenAGridLayoutPointerToBeOperational) -{ - TypeParam pair; - auto constexpr dim = pair.first(); - auto constexpr interp = pair.second(); - - using GridYee = GridLayout>; - - auto layout = std::make_unique(NDlayout::create()); - - // this->ohm.setLayout(layout.get()); - EXPECT_ANY_THROW( - this->ohm(this->n, this->V, this->P, this->B, this->J, - this->Enew)); // because the grid layout is not set (TODO issue #3392) -} - - std::vector read(std::string filename) { std::ifstream readFile(filename); @@ -338,6 +317,8 @@ std::vector read(std::string filename) TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) { + using Ohm_t = TestFixture::Ohm_t; + TypeParam pair; auto constexpr dim = pair.first(); auto constexpr interp = pair.second(); @@ -354,11 +335,9 @@ TYPED_TEST(OhmTest, ThatElectricFieldIsOkFromOhmsLaw) auto expected_ohmY = read(filenameY); auto expected_ohmZ = read(filenameZ); - using GridYee = GridLayout>; - auto layout = std::make_unique(NDlayout::create()); + auto const layout = NDlayout::create(); - this->ohm.setLayout(layout.get()); - this->ohm(this->n, this->V, this->P, this->B, this->J, this->Enew); + Ohm_t{this->ohm_info, layout}(this->n, this->V, this->P, this->B, this->J, this->Enew); auto const& [Exnew, Eynew, Eznew] = this->Enew(); From 9268d21189bb57468b8eb2e7601c236d4b98a0c8 Mon Sep 17 00:00:00 2001 From: deegan Date: Fri, 3 Apr 2026 14:37:14 +0200 Subject: [PATCH 3/3] level transformer --- .../hybrid_level_initializer.hpp | 14 +- .../resources_manager/resources_manager.hpp | 5 + src/amr/solvers/solver_field_evolvers.hpp | 134 ++++++++++++++++-- src/amr/solvers/solver_ppc.hpp | 49 ++++--- 4 files changed, 162 insertions(+), 40 deletions(-) diff --git a/src/amr/level_initializer/hybrid_level_initializer.hpp b/src/amr/level_initializer/hybrid_level_initializer.hpp index 2eb22896d..6dd0c0f9d 100644 --- a/src/amr/level_initializer/hybrid_level_initializer.hpp +++ b/src/amr/level_initializer/hybrid_level_initializer.hpp @@ -35,17 +35,15 @@ namespace solver static constexpr auto dimension = GridLayoutT::dimension; static constexpr auto interp_order = GridLayoutT::interp_order; - using Ampere_t = AmpereTransformer; - using Ohm_t = OhmTransformer; - - Ohm_t ohm_; - Ampere_t ampere_; + using Ampere_t = AmpereLevelTransformer; + using Ohm_t = OhmLevelTransformer; + core::OhmInfo ohm_info; inline bool isRootLevel(int const levelNumber) const { return levelNumber == 0; } public: explicit HybridLevelInitializer(PHARE::initializer::PHAREDict const& dict) - : ohm_{dict["algo"]["ohm"]} + : ohm_info{core::OhmInfo::FROM(dict["algo"]["ohm"])} { } @@ -158,7 +156,7 @@ namespace solver auto& E = hybridModel.state.electromag.E; auto& J = hybridModel.state.J; - ampere_(level, hybridModel, B, J); + Ampere_t{level, hybridModel}(B, J); setTime(J); hybMessenger.fillCurrentGhosts(J, level, 0.); @@ -166,7 +164,7 @@ namespace solver for (auto& patch : rm.enumerate(level, electrons)) electrons.update(amr::layoutFromPatch(*patch)); - ohm_(level, hybridModel, B, J, E); + Ohm_t{ohm_info, level, hybridModel}(B, J, E, electrons); setTime(E); hybMessenger.fillElectricGhosts(E, level, 0.); } diff --git a/src/amr/resources_manager/resources_manager.hpp b/src/amr/resources_manager/resources_manager.hpp index f660a2c92..87489f71e 100644 --- a/src/amr/resources_manager/resources_manager.hpp +++ b/src/amr/resources_manager/resources_manager.hpp @@ -362,6 +362,11 @@ namespace amr looper->set(*patch); } + Patch(Patch const&) = delete; + Patch(Patch&&) = delete; + Patch& operator=(Patch const&) = delete; + Patch& operator=(Patch&&) = delete; + SAMRAI::hier::Patch const& operator*() { return *patch; } ~Patch() { looper->unset(); } diff --git a/src/amr/solvers/solver_field_evolvers.hpp b/src/amr/solvers/solver_field_evolvers.hpp index 380281588..983307b3c 100644 --- a/src/amr/solvers/solver_field_evolvers.hpp +++ b/src/amr/solvers/solver_field_evolvers.hpp @@ -12,11 +12,12 @@ namespace PHARE::solver { -template +template class FaradayTransformer { - using core_type = core::Faraday; - using level_t = AMR_Types::level_t; + using GridLayout = Model::gridlayout_type; + using level_t = Model::amr_types::level_t; + using core_type = core::Faraday; public: void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } @@ -33,28 +34,101 @@ class FaradayTransformer } }; -template + +template +class FaradayLevelTransformer +{ + using GridLayout = Model::gridlayout_type; + using level_t = Model::amr_types::level_t; + using core_type = core::Faraday; + +public: + explicit FaradayLevelTransformer(level_t& level, auto& model) + : level_{level} + , model_{model} + { + } + + void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } + + void operator()(auto& B, auto& E, auto& Bnew, auto& dt) + { + auto& rm = *model_.resourcesManager; + for (auto& patch : rm.enumerate(level_, B, E, Bnew)) + { + auto layout = amr::layoutFromPatch(*patch); + (*this)(layout, B, E, Bnew, dt); + } + } + + level_t& level_; + Model& model_; +}; +template +FaradayLevelTransformer(typename Model::amr_types::level_t&, Model&) + -> FaradayLevelTransformer; + +template class AmpereTransformer { - using core_type = core::Ampere; - using level_t = AMR_Types::level_t; + using GridLayout = Model::gridlayout_type; + using level_t = Model::amr_types::level_t; + using core_type = core::Ampere; public: void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } void operator()(level_t& level, auto& model, auto& B, auto& J) { - auto& rm = *model.resourcesManager; - auto& state = model.state; - for (auto& patch : rm.enumerate(level, state, B, J)) + auto& rm = *model.resourcesManager; + for (auto& patch : rm.enumerate(level, B, J)) + { + auto layout = amr::layoutFromPatch(*patch); + (*this)(layout, B, J); + } + } +}; + + + + +template +class AmpereLevelTransformer +{ + using GridLayout = Model::gridlayout_type; + using level_t = Model::amr_types::level_t; + using core_type = core::Ampere; + +public: + explicit AmpereLevelTransformer(level_t& level, auto& model) + : level_{level} + , model_{model} + { + } + + void operator()(GridLayout& layout, auto&&... args) { core_type{layout}(args...); } + + void operator()(auto& B, auto& J) + { + auto& rm = *model_.resourcesManager; + for (auto& patch : rm.enumerate(level_, B, J)) { auto layout = amr::layoutFromPatch(*patch); (*this)(layout, B, J); } } + + level_t& level_; + Model& model_; }; +template +AmpereLevelTransformer(typename Model::amr_types::level_t&, Model&) + -> AmpereLevelTransformer; + + + template class OhmTransformer { @@ -93,6 +167,48 @@ class OhmTransformer }; +template +class OhmLevelTransformer +{ + using GridLayout = Model::gridlayout_type; + using level_t = Model::amr_types::level_t; + using info_type = core::OhmInfo; + using core_type = core::Ohm; + +public: + explicit OhmLevelTransformer(info_type const& info, level_t& level, Model& model) + : info_{info} + , level_{level} + , model_{model} + { + } + + void operator()(GridLayout& layout, auto&&... args) { core_type{info_, layout}(args...); } + + void operator()(auto& B, auto& J, auto& E, auto& electrons) + { + auto& rm = *model_.resourcesManager; + for (auto& patch : rm.enumerate(level_, electrons, B, J, E)) + { + auto layout = amr::layoutFromPatch(*patch); + auto& n = electrons.density(); + auto& Ve = electrons.velocity(); + auto& Pe = electrons.pressure(); + (*this)(layout, n, Ve, Pe, B, J, E); + } + } + + void operator()(auto& B, auto& E, auto& electrons) { (*this)(B, model_.state.J, E, electrons); } + + info_type info_; + level_t& level_; + Model& model_; +}; + +template +OhmLevelTransformer(core::OhmInfo, typename Model::amr_types::level_t&, Model&) + -> OhmLevelTransformer; + } // namespace PHARE::solver diff --git a/src/amr/solvers/solver_ppc.hpp b/src/amr/solvers/solver_ppc.hpp index 0b6b5cc25..6365b4ffc 100644 --- a/src/amr/solvers/solver_ppc.hpp +++ b/src/amr/solvers/solver_ppc.hpp @@ -1,14 +1,11 @@ #ifndef PHARE_SOLVER_PPC_HPP #define PHARE_SOLVER_PPC_HPP -#include "core/data/electrons/electrons.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep #include "core/numerics/ohm/ohm.hpp" #include "core/utilities/mpi_utils.hpp" -#include "core/numerics/ampere/ampere.hpp" #include "core/data/vecfield/vecfield.hpp" -#include "core/numerics/faraday/faraday.hpp" #include "core/numerics/ion_updater/ion_updater.hpp" #include "amr/solvers/solver.hpp" @@ -46,9 +43,9 @@ class SolverPPC : public ISolver using IMessenger = amr::IMessenger; using HybridMessenger = amr::HybridMessenger; - using Faraday_t = FaradayTransformer; - using Ampere_t = AmpereTransformer; - using Ohm_t = OhmTransformer; + using Faraday_t = FaradayLevelTransformer; + using Ampere_t = AmpereLevelTransformer; + using Ohm_t = OhmLevelTransformer; using IonUpdater_t = PHARE::core::IonUpdater; Electromag electromagPred_{"EMPred"}; @@ -58,10 +55,7 @@ class SolverPPC : public ISolver VecFieldT fluxSumE_{this->name() + "_fluxSumE", core::HybridQuantity::Vector::E}; std::unordered_map oldTime_; - Faraday_t faraday_; - Ampere_t ampere_; - Ohm_t ohm_; - + core::OhmInfo ohm_info; IonUpdater_t ionUpdater_; @@ -71,10 +65,9 @@ class SolverPPC : public ISolver using hierarchy_t = AMR_Types::hierarchy_t; - explicit SolverPPC(PHARE::initializer::PHAREDict const& dict) : ISolver{"PPC"} - , ohm_{dict["ohm"]} + , ohm_info{core::OhmInfo::FROM(dict["ohm"])} , ionUpdater_{dict["ion_updater"]} { @@ -323,7 +316,8 @@ void SolverPPC::reflux(IPhysicalModel_t& model, auto& Eavg = electromagAvg_.E; auto& B = hybridModel.state.electromag.B; auto dt = time - oldTime_[level.getLevelNumber()]; - faraday_(level, hybridModel, Bold_, Eavg, B, dt); + + Faraday_t{level, hybridModel}(Bold_, Eavg, B, dt); hybridMessenger.fillMagneticGhosts(B, level, time); } @@ -371,26 +365,29 @@ void SolverPPC::predictor1_(level_t& level, HybridModel& TimeSetter setTime{level, model, newTime}; + Faraday_t faraday{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.faraday"); auto dt = newTime - currentTime; - faraday_(level, model, model.state.electromag.B, model.state.electromag.E, - electromagPred_.B, dt); + faraday(model.state.electromag.B, model.state.electromag.E, electromagPred_.B, dt); setTime(electromagPred_.B); fromCoarser.fillMagneticGhosts(electromagPred_.B, level, newTime); } + Ampere_t ampere{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.ampere"); - ampere_(level, model, electromagPred_.B, model.state.J); + ampere(electromagPred_.B, model.state.J); setTime(model.state.J); fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } + + Ohm_t ohm{ohm_info, level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor1_.ohm"); update_electrons(level, model); - ohm_(level, model, electromagPred_.B, electromagPred_.E); + ohm(electromagPred_.B, electromagPred_.E, model.state.electrons); setTime(electromagPred_.E); } } @@ -405,25 +402,28 @@ void SolverPPC::predictor2_(level_t& level, HybridModel& TimeSetter setTime{level, model, newTime}; + Faraday_t faraday{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.faraday"); auto dt = newTime - currentTime; - faraday_(level, model, model.state.electromag.B, electromagAvg_.E, electromagPred_.B, dt); + faraday(model.state.electromag.B, electromagAvg_.E, electromagPred_.B, dt); setTime(electromagPred_.B); fromCoarser.fillMagneticGhosts(electromagPred_.B, level, newTime); } + Ampere_t ampere{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.ampere"); - ampere_(level, model, electromagPred_.B, model.state.J); + ampere(electromagPred_.B, model.state.J); setTime(model.state.J); fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } + Ohm_t ohm{ohm_info, level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::predictor2_.ohm"); update_electrons(level, model); - ohm_(level, model, electromagPred_.B, electromagPred_.E); + ohm(electromagPred_.B, electromagPred_.E, model.state.electrons); setTime(electromagPred_.E); } } @@ -442,26 +442,29 @@ void SolverPPC::corrector_(level_t& level, HybridModel& TimeSetter setTime{level, model, newTime}; auto& electromag = model.state.electromag; + Faraday_t faraday{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.faraday"); auto dt = newTime - currentTime; - faraday_(level, model, electromag.B, electromagAvg_.E, electromag.B, dt); + faraday(electromag.B, electromagAvg_.E, electromag.B, dt); setTime(model.state.electromag.B); fromCoarser.fillMagneticGhosts(model.state.electromag.B, level, newTime); } + Ampere_t ampere{level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.ampere"); - ampere_(level, model, electromag.B, model.state.J); + ampere(electromag.B, model.state.J); setTime(model.state.J); fromCoarser.fillCurrentGhosts(model.state.J, level, newTime); } + Ohm_t ohm{ohm_info, level, model}; { PHARE_LOG_SCOPE(1, "SolverPPC::corrector_.ohm"); update_electrons(level, model); - ohm_(level, model, electromag.B, electromag.E); + ohm(electromag.B, electromag.E, model.state.electrons); setTime(model.state.electromag.E); fromCoarser.fillElectricGhosts(model.state.electromag.E, level, newTime);