diff --git a/src/amr/solvers/solver_mhd.hpp b/src/amr/solvers/solver_mhd.hpp index ba2b78a8e..9a24f9369 100644 --- a/src/amr/solvers/solver_mhd.hpp +++ b/src/amr/solvers/solver_mhd.hpp @@ -16,7 +16,8 @@ #include "amr/solvers/solver_mhd_model_view.hpp" #include "core/data/grid/gridlayoutdefs.hpp" #include "core/data/vecfield/vecfield_component.hpp" -#include "core/numerics/mhd_reconstruction/mhd_reconstruction.hpp" +#include "core/mhd/mhd_quantities.hpp" +#include "core/numerics/godunov_fluxes/godunov_fluxes.hpp" namespace PHARE::solver { @@ -31,27 +32,41 @@ class SolverMHD : public ISolver using level_t = typename AMR_Types::level_t; using hierarchy_t = typename AMR_Types::hierarchy_t; - using field_type = typename MHDModel::field_type; - using GridLayout = typename MHDModel::gridlayout_type; + using FieldT = typename MHDModel::field_type; + using VecFieldT = typename MHDModel::vecfield_type; + using GridLayout = typename MHDModel::gridlayout_type; + using MHDQuantity = core::MHDQuantity; using IPhysicalModel_t = IPhysicalModel; using IMessenger = amr::IMessenger; using Direction = core::Direction; - using Reconstruction_t = core::Reconstruction; + using GodunovFluxes_t = typename MHDModelView::GodunovFluxes_t; + using Ampere_t = typename MHDModelView::Ampere_t; - std::vector uL_x; - std::vector uR_x; - std::vector uL_y; - std::vector uR_y; - std::vector uL_z; - std::vector uR_z; - Reconstruction_t reconstruct_; + FieldT rho_x{"rho_x", MHDQuantity::Scalar::ScalarFlux_x}; + VecFieldT rhoV_x{"V_x", MHDQuantity::Vector::VecFlux_x}; + VecFieldT B_x{"B_x", MHDQuantity::Vector::VecFlux_x}; + FieldT Etot_x{"rho_x", MHDQuantity::Scalar::ScalarFlux_x}; + + FieldT rho_y{"rho_y", MHDQuantity::Scalar::ScalarFlux_y}; + VecFieldT rhoV_y{"V_y", MHDQuantity::Vector::VecFlux_y}; + VecFieldT B_y{"B_y", MHDQuantity::Vector::VecFlux_y}; + FieldT Etot_y{"rho_y", MHDQuantity::Scalar::ScalarFlux_y}; + + FieldT rho_z{"rho_z", MHDQuantity::Scalar::ScalarFlux_z}; + VecFieldT rhoV_z{"V_z", MHDQuantity::Vector::VecFlux_z}; + VecFieldT B_z{"B_z", MHDQuantity::Vector::VecFlux_z}; + FieldT Etot_z{"rho_z", MHDQuantity::Scalar::ScalarFlux_z}; + + GodunovFluxes_t godunov_; + Ampere_t ampere_; public: - SolverMHD() + SolverMHD(PHARE::initializer::PHAREDict const& dict) : ISolver{"MHDSolver"} + , godunov_{dict["godunov"]} { } @@ -82,10 +97,7 @@ class SolverMHD : public ISolver } private: - void reconstruction_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, - double const currentTime, double const newTime); - - void riemann_solver_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, + void godunov_fluxes_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, double const currentTime, double const newTime); void FV_cycle_(level_t& level, ModelViews_t& views, Messenger& fromCoarser, @@ -123,9 +135,7 @@ void SolverMHD::advanceLevel( auto& fromCoarser = dynamic_cast(fromCoarserMessenger); auto level = hierarchy.getPatchLevel(levelNumber); - reconstruction_(*level, modelView, fromCoarser, currentTime, newTime); - - riemann_solver_(*level, modelView, fromCoarser, currentTime, newTime); + godunov_fluxes_(*level, modelView, fromCoarser, currentTime, newTime); FV_cycle_(*level, modelView, fromCoarser, currentTime, newTime); @@ -135,68 +145,34 @@ void SolverMHD::advanceLevel( } template -void SolverMHD::reconstruction_( +void SolverMHD::godunov_fluxes_( level_t& level, ModelViews_t& views, Messenger& fromCoarser, double const currentTime, double const newTime) { - PHARE_LOG_SCOPE(1, "SolverMHD::reconstruction_"); - - // Ampere - // centering - - auto Fields = std::forward_as_tuple( - views.super().rho, views.super().V(core::Component::X), views.super().V(core::Component::Y), - views.super().V(core::Component::Z), views.super().B_CT(core::Component::X), - views.super().B_CT(core::Component::Y), views.super().B_CT(core::Component::Z), - views.super().P); - - std::apply( - [&](auto&... field) { - if constexpr (dimension == 1) - { - [&](std::index_sequence) { - ((reconstruct_.template operator()(std::get(Fields), uL_x[I], - uR_x[I])), - ...); - }(std::make_index_sequence<8>{}); - } - if constexpr (dimension == 2) - { - [&](std::index_sequence) { - ((reconstruct_.template operator()(std::get(Fields), uL_x[I], - uR_x[I])), - ...); - ((reconstruct_.template operator()(std::get(Fields), uL_y[I], - uR_y[I])), - ...); - }(std::make_index_sequence<8>{}); - } - if constexpr (dimension == 3) - { - [&](std::index_sequence) { - ((reconstruct_.template operator()(std::get(Fields), uL_x[I], - uR_x[I])), - ...); - ((reconstruct_.template operator()(std::get(Fields), uL_y[I], - uR_y[I])), - ...); - ((reconstruct_.template operator()(std::get(Fields), uL_z[I], - uR_z[I])), - ...); - }(std::make_index_sequence<8>{}); - } - }, - Fields); -} + PHARE_LOG_SCOPE(1, "SolverMHD::godunov_fluxes_"); -template -void SolverMHD::riemann_solver_( - level_t& level, ModelViews_t& views, Messenger& fromCoarser, double const currentTime, - double const newTime) -{ - PHARE_LOG_SCOPE(1, "SolverMHD::riemann_solver_"); + ampere_(views.layouts, views.B_CT, views.J); + + if constexpr (dimension == 1) + { + godunov_(views.layouts, views.rho, views.V, views.B_CT, views.P, views.J, views.rho_x, + views.rhoV_x, views.B_x, views.Etot_x); + } + if constexpr (dimension == 2) + { + godunov_(views.layouts, views.rho, views.V, views.B_CT, views.P, views.J, views.rho_x, + views.rhoV_x, views.B_x, views.Etot_x, views.rho_y, views.rhoV_y, views.B_y, + views.Etot_y); + } + if constexpr (dimension == 3) + { + godunov_(views.layouts, views.rho, views.V, views.B_CT, views.P, views.J, views.rho_x, + views.rhoV_x, views.B_x, views.Etot_x, views.rho_y, views.rhoV_y, views.B_y, + views.Etot_y, views.rho_z, views.rhoV_z, views.B_z, views.Etot_z); + } } + template void SolverMHD::FV_cycle_(level_t& level, ModelViews_t& views, @@ -221,6 +197,7 @@ void SolverMHD::constrained_transp auto dt = newTime - currentTime; + // Fill ghost for B_RS // averaging B_RS(x, y, z) // faraday } diff --git a/src/amr/solvers/solver_mhd_model_view.hpp b/src/amr/solvers/solver_mhd_model_view.hpp index aea02c15c..d066d5407 100644 --- a/src/amr/solvers/solver_mhd_model_view.hpp +++ b/src/amr/solvers/solver_mhd_model_view.hpp @@ -2,12 +2,43 @@ #define PHARE_SOLVER_SOLVER_MHD_MODEL_VIEW_HPP #include "amr/solvers/solver.hpp" +#include "amr/solvers/solver_ppc_model_view.hpp" +#include "core/numerics/godunov_fluxes/godunov_fluxes.hpp" -namespace PHARE::solver { +namespace PHARE::solver +{ +template +class GodunovFluxesTransformer +{ + using core_type = PHARE::core::GodunovFluxes; -template -class MHDModelView : public ISolverModelView {}; +public: + template + void operator()(Layout const& layouts, Field const& rho, VecField const& V, + VecField const& B_CT, Field const& P, VecField const& J, Fluxes&... fluxes) + { + assert_equal_sizes(rho, V, B_CT, P, J, fluxes...); + for (std::size_t i = 0; i < layouts.size(); ++i) + { + auto _ = core::SetLayout(layouts[i], godunov_); + godunov_(*rho[i], *V[i], *B_CT[i], *P[i], *J[i], *fluxes[i]...); + } + } -}; // namespace PHARE::solver + core_type godunov_; +}; -#endif // PHARE_SOLVER_SOLVER_MHD_MODEL_VIEW_HPP + +template +class MHDModelView : public ISolverModelView +{ +public: + using MHDModel_t = MHDModel_; + using GridLayout = typename MHDModel_t::gridlayout_type; + using GodunovFluxes_t = GodunovFluxesTransformer; + using Ampere_t = AmpereTransformer; +}; + +}; // namespace PHARE::solver + +#endif // PHARE_SOLVER_SOLVER_MHD_MODEL_VIEW_HPP diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index ac9929742..1d671c8d7 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -287,7 +287,7 @@ namespace core return physicalStartIndexTable_[icentering][iDir]; } - NO_DISCARD std::uint32_t physicalStartIndex(Quantity::Scalar const& quantity, + NO_DISCARD std::uint32_t physicalStartIndex(typename Quantity::Scalar const& quantity, Direction direction) const { std::uint32_t iQty = static_cast(quantity); @@ -322,7 +322,7 @@ namespace core return physicalEndIndexTable_[icentering][iDir]; } - NO_DISCARD std::uint32_t physicalEndIndex(Quantity::Scalar const& quantity, + NO_DISCARD std::uint32_t physicalEndIndex(typename Quantity::Scalar const& quantity, Direction direction) const { std::uint32_t iQty = static_cast(quantity); @@ -357,7 +357,8 @@ namespace core return 0; } - NO_DISCARD std::uint32_t ghostStartIndex([[maybe_unused]] Quantity::Scalar const& quantity, + NO_DISCARD std::uint32_t ghostStartIndex([[maybe_unused]] + typename Quantity::Scalar const& quantity, [[maybe_unused]] Direction direction) const { // ghostStartIndex is always the first node @@ -389,7 +390,7 @@ namespace core return ghostEndIndexTable_[iCentering][iDir]; } - NO_DISCARD std::uint32_t ghostEndIndex(Quantity::Scalar const& quantity, + NO_DISCARD std::uint32_t ghostEndIndex(typename Quantity::Scalar const& quantity, Direction direction) const { std::uint32_t iQty = static_cast(quantity); @@ -809,7 +810,7 @@ namespace core * @brief returns the centering of a scalar hybrid quantity in each directions */ NO_DISCARD constexpr static std::array - centering(Quantity::Scalar quantity) + centering(typename Quantity::Scalar quantity) { return GridLayoutImpl::centering(quantity); } @@ -818,7 +819,7 @@ namespace core * @brief returns the centering of a vector hybrid quantity in each directions */ NO_DISCARD constexpr static std::array, 3> - centering(Quantity::Vector quantity) + centering(typename Quantity::Vector quantity) { return GridLayoutImpl::centering(quantity); } @@ -828,7 +829,8 @@ namespace core * @return An std::array object, containing the size to which allocate * arrays of an Quantity::Quantity 'qty' in every directions. */ - NO_DISCARD std::array allocSize(Quantity::Scalar qty) const + NO_DISCARD std::array + allocSize(typename Quantity::Scalar qty) const { std::uint32_t iQty = static_cast(qty); @@ -851,8 +853,8 @@ namespace core * @brief allocSizeDerived returns the shape of the array to be allocated to store * the derivative of a given quantity in a given direction. */ - NO_DISCARD std::array allocSizeDerived(Quantity::Scalar qty, - Direction dir) const + NO_DISCARD std::array + allocSizeDerived(typename Quantity::Scalar qty, Direction dir) const { std::uint32_t iDerivedDir = static_cast(dir); std::uint32_t iQty = static_cast(qty); @@ -893,7 +895,7 @@ namespace core * ghost nodes */ NO_DISCARD std::array - nbrPhysicalNodes(Quantity::Scalar hybQty) const + nbrPhysicalNodes(typename Quantity::Scalar hybQty) const { std::array centerings; @@ -912,7 +914,7 @@ namespace core * primal and primal becomes dual. quantityCentering is used to know if the * Quantity::Quantity 'qty' is primal or dual in the Direction 'dir' */ - NO_DISCARD QtyCentering derivedCentering(Quantity::Scalar qty, Direction dir) const + NO_DISCARD QtyCentering derivedCentering(typename Quantity::Scalar qty, Direction dir) const { std::uint32_t iField = static_cast(qty); std::uint32_t idir = static_cast(dir); diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp index bfcb06f30..208cee720 100644 --- a/src/core/data/grid/gridlayoutdefs.hpp +++ b/src/core/data/grid/gridlayoutdefs.hpp @@ -141,26 +141,33 @@ namespace core static constexpr std::uint32_t iJy = static_cast(MHDQuantity::Scalar::Jy); static constexpr std::uint32_t iJz = static_cast(MHDQuantity::Scalar::Jz); - static constexpr std::uint32_t iBx_RSx - = static_cast(MHDQuantity::Scalar::Bx_RSx); - static constexpr std::uint32_t iBy_RSx - = static_cast(MHDQuantity::Scalar::By_RSx); - static constexpr std::uint32_t iBz_RSx - = static_cast(MHDQuantity::Scalar::Bz_RSx); - - static constexpr std::uint32_t iBx_RSy - = static_cast(MHDQuantity::Scalar::Bx_RSy); - static constexpr std::uint32_t iBy_RSy - = static_cast(MHDQuantity::Scalar::By_RSy); - static constexpr std::uint32_t iBz_RSy - = static_cast(MHDQuantity::Scalar::Bz_RSy); - - static constexpr std::uint32_t iBx_RSz - = static_cast(MHDQuantity::Scalar::Bx_RSz); - static constexpr std::uint32_t iBy_RSz - = static_cast(MHDQuantity::Scalar::By_RSz); - static constexpr std::uint32_t iBz_RSz - = static_cast(MHDQuantity::Scalar::Bz_RSz); + static constexpr std::uint32_t iScalarFlux_x + = static_cast(MHDQuantity::Scalar::ScalarFlux_x); + static constexpr std::uint32_t iScalarFlux_y + = static_cast(MHDQuantity::Scalar::ScalarFlux_y); + static constexpr std::uint32_t iScalarFlux_z + = static_cast(MHDQuantity::Scalar::ScalarFlux_z); + + static constexpr std::uint32_t iVecFluxX_x + = static_cast(MHDQuantity::Scalar::VecFluxX_x); + static constexpr std::uint32_t iVecFluxY_x + = static_cast(MHDQuantity::Scalar::VecFluxY_x); + static constexpr std::uint32_t iVecFluxZ_x + = static_cast(MHDQuantity::Scalar::VecFluxZ_x); + + static constexpr std::uint32_t iVecFluxX_y + = static_cast(MHDQuantity::Scalar::VecFluxX_y); + static constexpr std::uint32_t iVecFluxY_y + = static_cast(MHDQuantity::Scalar::VecFluxY_y); + static constexpr std::uint32_t iVecFluxZ_y + = static_cast(MHDQuantity::Scalar::VecFluxZ_y); + + static constexpr std::uint32_t iVecFluxX_z + = static_cast(MHDQuantity::Scalar::VecFluxX_z); + static constexpr std::uint32_t iVecFluxY_z + = static_cast(MHDQuantity::Scalar::VecFluxY_z); + static constexpr std::uint32_t iVecFluxZ_z + = static_cast(MHDQuantity::Scalar::VecFluxZ_z); }; } // namespace core } // namespace PHARE diff --git a/src/core/data/grid/gridlayoutimplyee_mhd.hpp b/src/core/data/grid/gridlayoutimplyee_mhd.hpp index 21ccf9f37..2a5c80276 100644 --- a/src/core/data/grid/gridlayoutimplyee_mhd.hpp +++ b/src/core/data/grid/gridlayoutimplyee_mhd.hpp @@ -84,34 +84,46 @@ namespace core const std::array Jz = {{data.primal, data.primal, data.dual}}; - const std::array Bx_RSx + const std::array ScalarFlux_x = {{data.primal, data.dual, data.dual}}; - const std::array By_RSx + + const std::array ScalarFlux_y + = {{data.dual, data.primal, data.dual}}; + + const std::array ScalarFlux_z + = {{data.dual, data.dual, data.primal}}; + + const std::array VecFluxX_x = {{data.primal, data.dual, data.dual}}; - const std::array Bz_RSx + const std::array VecFluxY_x + = {{data.primal, data.dual, data.dual}}; + const std::array VecFluxZ_x = {{data.primal, data.dual, data.dual}}; - const std::array Bx_RSy + const std::array VecFluxX_y = {{data.dual, data.primal, data.dual}}; - const std::array By_RSy + const std::array VecFluxY_y = {{data.dual, data.primal, data.dual}}; - const std::array Bz_RSy + const std::array VecFluxZ_y = {{data.dual, data.primal, data.dual}}; - const std::array Bx_RSz + const std::array VecFluxX_z = {{data.dual, data.dual, data.primal}}; - const std::array By_RSz + const std::array VecFluxY_z = {{data.dual, data.dual, data.primal}}; - const std::array Bz_RSz + const std::array VecFluxZ_z = {{data.dual, data.dual, data.primal}}; const std::array, static_cast(MHDQuantity::Scalar::count)> - _QtyCentering{Rho, Vx, Vy, Vz, P, rhoVx, rhoVy, rhoVz, - Bx_FV, By_FV, Bz_FV, Etot, Bx_CT, By_CT, Bz_CT, Ex, - Ey, Ez, Jx, Jy, Jz, Bx_RSx, By_RSx, Bz_RSx, - Bx_RSy, By_RSy, Bz_RSy, Bx_RSz, By_RSz, Bz_RSz}; + _QtyCentering{Rho, Vx, Vy, Vz, P, + rhoVx, rhoVy, rhoVz, Bx_FV, By_FV, + Bz_FV, Etot, Bx_CT, By_CT, Bz_CT, + Ex, Ey, Ez, Jx, Jy, + Jz, ScalarFlux_x, ScalarFlux_y, ScalarFlux_z, VecFluxX_x, + VecFluxY_x, VecFluxZ_x, VecFluxX_y, VecFluxY_y, VecFluxZ_y, + VecFluxX_z, VecFluxY_z, VecFluxZ_z}; return _QtyCentering; } @@ -175,24 +187,14 @@ namespace core return {{_QtyCentering_[gridData_.iJy][gridData_.idirX]}}; case MHDQuantity::Scalar::Jz: return {{_QtyCentering_[gridData_.iJz][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bx_RSx: - return {{_QtyCentering_[gridData_.iBx_RSx][gridData_.idirX]}}; - case MHDQuantity::Scalar::By_RSx: - return {{_QtyCentering_[gridData_.iBy_RSx][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bz_RSx: - return {{_QtyCentering_[gridData_.iBz_RSx][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bx_RSy: - return {{_QtyCentering_[gridData_.iBx_RSy][gridData_.idirX]}}; - case MHDQuantity::Scalar::By_RSy: - return {{_QtyCentering_[gridData_.iBy_RSy][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bz_RSy: - return {{_QtyCentering_[gridData_.iBz_RSy][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bx_RSz: - return {{_QtyCentering_[gridData_.iBx_RSz][gridData_.idirX]}}; - case MHDQuantity::Scalar::By_RSz: - return {{_QtyCentering_[gridData_.iBy_RSz][gridData_.idirX]}}; - case MHDQuantity::Scalar::Bz_RSz: - return {{_QtyCentering_[gridData_.iBz_RSz][gridData_.idirX]}}; + case MHDQuantity::Scalar::ScalarFlux_x: + return {{_QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirX]}}; + case MHDQuantity::Scalar::VecFluxX_x: + return {{_QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirX]}}; + case MHDQuantity::Scalar::VecFluxY_x: + return {{_QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirX]}}; + case MHDQuantity::Scalar::VecFluxZ_x: + return {{_QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirX]}}; default: throw std::runtime_error("Wrong MHDQuantity"); } } @@ -261,33 +263,30 @@ namespace core case MHDQuantity::Scalar::Jz: return {{_QtyCentering_[gridData_.iJz][gridData_.idirX], _QtyCentering_[gridData_.iJz][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bx_RSx: - return {{_QtyCentering_[gridData_.iBx_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSx][gridData_.idirY]}}; - case MHDQuantity::Scalar::By_RSx: - return {{_QtyCentering_[gridData_.iBy_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSx][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bz_RSx: - return {{_QtyCentering_[gridData_.iBz_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSx][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bx_RSy: - return {{_QtyCentering_[gridData_.iBx_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSy][gridData_.idirY]}}; - case MHDQuantity::Scalar::By_RSy: - return {{_QtyCentering_[gridData_.iBy_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSy][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bz_RSy: - return {{_QtyCentering_[gridData_.iBz_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSy][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bx_RSz: - return {{_QtyCentering_[gridData_.iBx_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSz][gridData_.idirY]}}; - case MHDQuantity::Scalar::By_RSz: - return {{_QtyCentering_[gridData_.iBy_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSz][gridData_.idirY]}}; - case MHDQuantity::Scalar::Bz_RSz: - return {{_QtyCentering_[gridData_.iBz_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSz][gridData_.idirY]}}; + case MHDQuantity::Scalar::ScalarFlux_x: + return {{_QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirX], + _QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirY]}}; + case MHDQuantity::Scalar::ScalarFlux_y: + return {{_QtyCentering_[gridData_.iScalarFlux_y][gridData_.idirX], + _QtyCentering_[gridData_.iScalarFlux_y][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxX_x: + return {{_QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxY_x: + return {{_QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxZ_x: + return {{_QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxX_y: + return {{_QtyCentering_[gridData_.iVecFluxX_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxX_y][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxY_y: + return {{_QtyCentering_[gridData_.iVecFluxY_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxY_y][gridData_.idirY]}}; + case MHDQuantity::Scalar::VecFluxZ_y: + return {{_QtyCentering_[gridData_.iVecFluxZ_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxZ_y][gridData_.idirY]}}; default: throw std::runtime_error("Wrong MHDQuantity"); } } @@ -376,42 +375,54 @@ namespace core return {{_QtyCentering_[gridData_.iJz][gridData_.idirX], _QtyCentering_[gridData_.iJz][gridData_.idirY], _QtyCentering_[gridData_.iJz][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bx_RSx: - return {{_QtyCentering_[gridData_.iBx_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSx][gridData_.idirY], - _QtyCentering_[gridData_.iBx_RSx][gridData_.idirZ]}}; - case MHDQuantity::Scalar::By_RSx: - return {{_QtyCentering_[gridData_.iBy_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSx][gridData_.idirY], - _QtyCentering_[gridData_.iBy_RSx][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bz_RSx: - return {{_QtyCentering_[gridData_.iBz_RSx][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSx][gridData_.idirY], - _QtyCentering_[gridData_.iBz_RSx][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bx_RSy: - return {{_QtyCentering_[gridData_.iBx_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSy][gridData_.idirY], - _QtyCentering_[gridData_.iBx_RSy][gridData_.idirZ]}}; - case MHDQuantity::Scalar::By_RSy: - return {{_QtyCentering_[gridData_.iBy_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSy][gridData_.idirY], - _QtyCentering_[gridData_.iBy_RSy][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bz_RSy: - return {{_QtyCentering_[gridData_.iBz_RSy][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSy][gridData_.idirY], - _QtyCentering_[gridData_.iBz_RSy][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bx_RSz: - return {{_QtyCentering_[gridData_.iBx_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBx_RSz][gridData_.idirY], - _QtyCentering_[gridData_.iBx_RSz][gridData_.idirZ]}}; - case MHDQuantity::Scalar::By_RSz: - return {{_QtyCentering_[gridData_.iBy_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBy_RSz][gridData_.idirY], - _QtyCentering_[gridData_.iBy_RSz][gridData_.idirZ]}}; - case MHDQuantity::Scalar::Bz_RSz: - return {{_QtyCentering_[gridData_.iBz_RSz][gridData_.idirX], - _QtyCentering_[gridData_.iBz_RSz][gridData_.idirY], - _QtyCentering_[gridData_.iBz_RSz][gridData_.idirZ]}}; + case MHDQuantity::Scalar::ScalarFlux_x: + return {{_QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirX], + _QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirY], + _QtyCentering_[gridData_.iScalarFlux_x][gridData_.idirZ]}}; + case MHDQuantity::Scalar::ScalarFlux_y: + return {{_QtyCentering_[gridData_.iScalarFlux_y][gridData_.idirX], + _QtyCentering_[gridData_.iScalarFlux_y][gridData_.idirY], + _QtyCentering_[gridData_.iScalarFlux_y][gridData_.idirZ]}}; + case MHDQuantity::Scalar::ScalarFlux_z: + return {{_QtyCentering_[gridData_.iScalarFlux_z][gridData_.idirX], + _QtyCentering_[gridData_.iScalarFlux_z][gridData_.idirY], + _QtyCentering_[gridData_.iScalarFlux_z][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxX_x: + return {{_QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxX_x][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxY_x: + return {{_QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxY_x][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxZ_x: + return {{_QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxZ_x][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxX_y: + return {{_QtyCentering_[gridData_.iVecFluxX_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxX_y][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxX_y][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxY_y: + return {{_QtyCentering_[gridData_.iVecFluxY_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxY_y][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxY_y][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxZ_y: + return {{_QtyCentering_[gridData_.iVecFluxZ_y][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxZ_y][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxZ_y][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxX_z: + return {{_QtyCentering_[gridData_.iVecFluxX_z][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxX_z][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxX_z][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxY_z: + return {{_QtyCentering_[gridData_.iVecFluxY_z][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxY_z][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxY_z][gridData_.idirZ]}}; + case MHDQuantity::Scalar::VecFluxZ_z: + return {{_QtyCentering_[gridData_.iVecFluxZ_z][gridData_.idirX], + _QtyCentering_[gridData_.iVecFluxZ_z][gridData_.idirY], + _QtyCentering_[gridData_.iVecFluxZ_z][gridData_.idirZ]}}; default: throw std::runtime_error("Wrong MHDQuantity"); } } @@ -449,20 +460,20 @@ namespace core return {{centering(MHDQuantity::Scalar::Jx), centering(MHDQuantity::Scalar::Jy), centering(MHDQuantity::Scalar::Jz)}}; - case MHDQuantity::Vector::B_RSx: - return {{centering(MHDQuantity::Scalar::Bx_RSx), - centering(MHDQuantity::Scalar::By_RSx), - centering(MHDQuantity::Scalar::Bz_RSx)}}; + case MHDQuantity::Vector::VecFlux_x: + return {{centering(MHDQuantity::Scalar::VecFluxX_x), + centering(MHDQuantity::Scalar::VecFluxY_x), + centering(MHDQuantity::Scalar::VecFluxZ_x)}}; - case MHDQuantity::Vector::B_RSy: - return {{centering(MHDQuantity::Scalar::Bx_RSy), - centering(MHDQuantity::Scalar::By_RSy), - centering(MHDQuantity::Scalar::Bz_RSy)}}; + case MHDQuantity::Vector::VecFlux_y: + return {{centering(MHDQuantity::Scalar::VecFluxX_y), + centering(MHDQuantity::Scalar::VecFluxY_y), + centering(MHDQuantity::Scalar::VecFluxZ_y)}}; - case MHDQuantity::Vector::B_RSz: - return {{centering(MHDQuantity::Scalar::Bx_RSz), - centering(MHDQuantity::Scalar::By_RSz), - centering(MHDQuantity::Scalar::Bz_RSz)}}; + case MHDQuantity::Vector::VecFlux_z: + return {{centering(MHDQuantity::Scalar::VecFluxX_z), + centering(MHDQuantity::Scalar::VecFluxY_z), + centering(MHDQuantity::Scalar::VecFluxZ_z)}}; default: throw std::runtime_error("Wrong MHDQuantity"); } @@ -471,6 +482,102 @@ namespace core NO_DISCARD auto static constexpr dualToPrimal() { return -1; } NO_DISCARD auto static constexpr primalToDual() { return 1; } + + /* + NO_DISCARD auto static constexpr quantitiesToFaceX() + { + // The mhd quantities in FV are Ddd + // the X face is Pdd + // operation is thus Ddd to Pdd + // shift only in the X direction + + auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + constexpr WeightPoint P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{iShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.5}; + constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; + return std::array, 2>{P1, P2}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; + constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; + return std::array, 2>{P1, P2}; + } + } + + NO_DISCARD auto static constexpr quantitiesToFaceY() + { + // The mhd quantities in FV are Ddd + // the Y face is Dpd + // operation is thus Ddd to Dpd + // shift only in the Y direction + + [[maybe_unused]] auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + // since the linear combination is in the Y direction + // in 1D the quantities are already on the Y face so return 1 point with no + shift + // with coef 1. + constexpr WeightPoint P1{Point{0}, 1}; + return std::array, 1>{P1}; + } + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.5}; + constexpr WeightPoint P2{Point{0, iShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; + constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; + return std::array, 2>{P1, P2}; + } + } + + NO_DISCARD auto static constexpr quantitiesToFaceZ() + { + // The mhd quantities in FV are Ddd + // the Z face is Ddp + // operation is thus Ddd to Ddp + // shift only in the Z direction + + [[maybe_unused]] auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + // since the linear combination is in the Z direction + // in 1D or 2D the quantities are already on the Z face so return 1 point with + no + // shift with coef 1. + constexpr WeightPoint P1{Point{0}, 1.}; + return std::array, 1>{P1}; + } + else if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 1.}; + return std::array, 1>{P1}; + } + else if constexpr (dimension == 3) + { + // in 3D we need two points, the second with a primalToDual shift along Z + constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; + constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; + return std::array, 2>{P1, P2}; + } + } + */ + }; // namespace core } // namespace core diff --git a/src/core/mhd/mhd_quantities.hpp b/src/core/mhd/mhd_quantities.hpp index e57231c45..9cb301560 100644 --- a/src/core/mhd/mhd_quantities.hpp +++ b/src/core/mhd/mhd_quantities.hpp @@ -37,21 +37,22 @@ class MHDQuantity Jy, Jz, - Bx_RSx, - By_RSx, - Bz_RSx, - - Bx_RSy, - By_RSy, - Bz_RSy, - - Bx_RSz, - By_RSz, - Bz_RSz, + ScalarFlux_x, + ScalarFlux_y, + ScalarFlux_z, + VecFluxX_x, + VecFluxY_x, + VecFluxZ_x, + VecFluxX_y, + VecFluxY_y, + VecFluxZ_y, + VecFluxX_z, + VecFluxY_z, + VecFluxZ_z, count }; - enum class Vector { V, B_FV, rhoV, B_CT, E, J, B_RSx, B_RSy, B_RSz }; + enum class Vector { V, B_FV, rhoV, B_CT, E, J, VecFlux_x, VecFlux_y, VecFlux_z }; enum class Tensor { count }; template> @@ -65,9 +66,9 @@ class MHDQuantity NO_DISCARD static constexpr auto E() { return componentsQuantities(Vector::E); } NO_DISCARD static constexpr auto J() { return componentsQuantities(Vector::J); } - NO_DISCARD static constexpr auto B_RSx() { return componentsQuantities(Vector::B_RSx); } - NO_DISCARD static constexpr auto B_RSy() { return componentsQuantities(Vector::B_RSy); } - NO_DISCARD static constexpr auto B_RSz() { return componentsQuantities(Vector::B_RSz); } + NO_DISCARD static constexpr auto VecFlux_x() { return componentsQuantities(Vector::VecFlux_x); } + NO_DISCARD static constexpr auto VecFlux_y() { return componentsQuantities(Vector::VecFlux_y); } + NO_DISCARD static constexpr auto VecFlux_z() { return componentsQuantities(Vector::VecFlux_z); } NO_DISCARD static constexpr std::array componentsQuantities(Vector qty) { @@ -91,14 +92,14 @@ class MHDQuantity return {{Scalar::Jx, Scalar::Jy, Scalar::Jz}}; - if (qty == Vector::B_RSx) - return {{Scalar::Bx_RSx, Scalar::By_RSx, Scalar::Bz_RSx}}; + if (qty == Vector::VecFlux_x) + return {{Scalar::VecFluxX_x, Scalar::VecFluxY_x, Scalar::VecFluxZ_x}}; - if (qty == Vector::B_RSy) - return {{Scalar::Bx_RSy, Scalar::By_RSy, Scalar::Bz_RSy}}; + if (qty == Vector::VecFlux_y) + return {{Scalar::VecFluxX_y, Scalar::VecFluxY_y, Scalar::VecFluxZ_y}}; - if (qty == Vector::B_RSz) - return {{Scalar::Bx_RSz, Scalar::By_RSz, Scalar::Bz_RSz}}; + if (qty == Vector::VecFlux_z) + return {{Scalar::VecFluxX_z, Scalar::VecFluxY_z, Scalar::VecFluxZ_z}}; throw std::runtime_error("Error - invalid Vector"); } diff --git a/src/core/models/mhd_state.hpp b/src/core/models/mhd_state.hpp index 0b896d1a8..6f01ce4e1 100644 --- a/src/core/models/mhd_state.hpp +++ b/src/core/models/mhd_state.hpp @@ -29,27 +29,24 @@ namespace core { return rho.isUsable() and V.isUsable() and B_FV.isUsable() and P.isUsable() and rhoV.isUsable() and Etot.isUsable() and B_CT.isUsable() and J.isUsable() - and E.isUsable() and B_RSx.isUsable() and B_RSy.isUsable() and B_RSz.isUsable(); + and E.isUsable(); } NO_DISCARD bool isSettable() const { return rho.isSettable() and V.isSettable() and B_FV.isSettable() and P.isSettable() and rhoV.isSettable() and Etot.isSettable() and B_CT.isSettable() - and J.isSettable() and E.isSettable() and B_RSx.isSettable() - and B_RSy.isSettable() and B_RSz.isSettable(); + and J.isSettable() and E.isSettable(); } NO_DISCARD auto getCompileTimeResourcesViewList() const { - return std::forward_as_tuple(rho, V, B_FV, P, rhoV, Etot, B_CT, J, E, B_RSx, B_RSy, - B_RSz); + return std::forward_as_tuple(rho, V, B_FV, P, rhoV, Etot, B_CT, J, E); } NO_DISCARD auto getCompileTimeResourcesViewList() { - return std::forward_as_tuple(rho, V, B_FV, P, rhoV, Etot, B_CT, J, E, B_RSx, B_RSy, - B_RSz); + return std::forward_as_tuple(rho, V, B_FV, P, rhoV, Etot, B_CT, J, E); } //------------------------------------------------------------------------- @@ -72,13 +69,6 @@ namespace core , J{"J", MHDQuantity::Vector::J} , - B_RSx{"B_RS", MHDQuantity::Vector::B_RSx} - , B_RSy{"B_RS", MHDQuantity::Vector::B_RSy} - , B_RSz{"B_RS", MHDQuantity::Vector::B_RSz} - , - - - rhoinit_{ dict["density"]["initializer"].template to>()} , Vinit_{dict["velocity"]["initializer"]} @@ -109,11 +99,6 @@ namespace core VecFieldT E; VecFieldT J; - VecFieldT B_RSx; - VecFieldT B_RSy; - VecFieldT B_RSz; - - private: initializer::InitFunction rhoinit_; VecFieldInitializer Vinit_; diff --git a/src/core/numerics/godunov_fluxes/godunov_fluxes.hpp b/src/core/numerics/godunov_fluxes/godunov_fluxes.hpp new file mode 100644 index 000000000..467d3982e --- /dev/null +++ b/src/core/numerics/godunov_fluxes/godunov_fluxes.hpp @@ -0,0 +1,594 @@ +#ifndef PHARE_CORE_NUMERICS_GODUNOV_FLUXES_HPP +#define PHARE_CORE_NUMERICS_GODUNOV_FLUXES_HPP + +#include "core/data/grid/gridlayout.hpp" +#include "core/data/grid/gridlayout_utils.hpp" +#include "core/data/grid/gridlayoutdefs.hpp" +#include "core/data/vecfield/vecfield_component.hpp" +#include "core/utilities/index/index.hpp" +#include "initializer/data_provider.hpp" +#include +#include + +namespace PHARE::core +{ +template +class GodunovFluxes : public LayoutHolder +{ + constexpr static auto dimension = GridLayout::dimension; + using LayoutHolder::layout_; + +public: + GodunovFluxes(PHARE::initializer::PHAREDict const& dict) + : eta_{dict["resistivity"].template to()} + , nu_{dict["hyper_resistivity"].template to()} + , gamma_{dict["heat_capacity_ratio"].template to()} + { + } + + template + void operator()(Field const& rho, VecField const& V, VecField const& B_CT, Field const& P, + VecField const& J, Fluxes&... fluxes) + { + if (!this->hasLayout()) + throw std::runtime_error("Error - Reconstruction - GridLayout not set, cannot proceed " + "to reconstruction"); + + if constexpr (dimension == 1) + { + auto&& flux_tuple = std::forward_as_tuple(fluxes...); + + auto& rho_x = std::get<0>(flux_tuple); + auto& rhoV_x = std::get<1>(flux_tuple); + auto& B_x = std::get<2>(flux_tuple); + auto& Etot_x = std::get<3>(flux_tuple); + + layout_->evalOnBox(rho_x, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_x, rhoV_x, B_x, + Etot_x, {args...}); + }); + } + if constexpr (dimension == 2) + { + auto&& flux_tuple = std::forward_as_tuple(fluxes...); + + auto& rho_x = std::get<0>(flux_tuple); + auto& rhoV_x = std::get<1>(flux_tuple); + auto& B_x = std::get<2>(flux_tuple); + auto& Etot_x = std::get<3>(flux_tuple); + + auto& rho_y = std::get<4>(flux_tuple); + auto& rhoV_y = std::get<5>(flux_tuple); + auto& B_y = std::get<6>(flux_tuple); + auto& Etot_y = std::get<7>(flux_tuple); + + layout_->evalOnBox(rho_x, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_x, rhoV_x, B_x, + Etot_x, {args...}); + }); + + layout_->evalOnBox(rho_y, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_y, rhoV_y, B_y, + Etot_y, {args...}); + }); + } + if constexpr (dimension == 3) + { + auto&& flux_tuple = std::forward_as_tuple(fluxes...); + + auto& rho_x = std::get<0>(flux_tuple); + auto& rhoV_x = std::get<1>(flux_tuple); + auto& B_x = std::get<2>(flux_tuple); + auto& Etot_x = std::get<3>(flux_tuple); + + auto& rho_y = std::get<4>(flux_tuple); + auto& rhoV_y = std::get<5>(flux_tuple); + auto& B_y = std::get<6>(flux_tuple); + auto& Etot_y = std::get<7>(flux_tuple); + + auto& rho_z = std::get<8>(flux_tuple); + auto& rhoV_z = std::get<9>(flux_tuple); + auto& B_z = std::get<10>(flux_tuple); + auto& Etot_z = std::get<11>(flux_tuple); + + layout_->evalOnBox(rho_x, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_x, rhoV_x, B_x, + Etot_x, {args...}); + }); + + layout_->evalOnBox(rho_y, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_y, rhoV_y, B_y, + Etot_y, {args...}); + }); + + layout_->evalOnBox(rho_z, [&](auto&... args) mutable { + this->template godunov_fluxes_(rho, V, B_CT, P, J, rho_z, rhoV_z, B_z, + Etot_z, {args...}); + }); + } + } + +private: + double const gamma_; + double const eta_; + double const nu_; + + template + void godunov_fluxes_(Field const& rho, VecField const& V, VecField const& B_CT, Field const& P, + VecField const& J, Field& rho_flux, VecField& rhoV_flux, VecField& B_flux, + Field& Etot_flux, MeshIndex index) const + { + auto const& [Vx, Vy, Vz] = V(); + auto const& [Bx, By, Bz] = B_CT(); + auto const& [Jx, Jy, Jz] = J(); + + auto const& Quantities = std::forward_as_tuple(rho, Vx, Vy, Vz, Bx, By, Bz, P, Jx, Jy, Jz); + + // Left and Right state reconstructions + auto [rhoL, VxL, VyL, VzL, BxL, ByL, BzL, PL, JxL, JyL, JzL] = std::apply( + [&](auto const&... fields) { + return std::make_tuple(reconstruct_uL_(fields, index)...); + }, + Quantities); + + auto [rhoR, VxR, VyR, VzR, BxR, ByR, BzR, PR, JxR, JyR, JzR] = std::apply( + [&](auto const&... fields) { + return std::make_tuple(reconstruct_uR_(fields, index)...); + }, + Quantities); + + rhoR = std::abs(rhoR); + + // Compute ideal flux vector for Left and Right states + auto [F_rhoL, F_rhoVxL, F_rhoVyL, F_rhoVzL, F_BxL, F_ByL, F_BzL, F_EtotL] + = ideal_flux_vector_(rhoL, VxL, VyL, VzL, BxL, ByL, BzL, PL); + + auto [F_rhoR, F_rhoVxR, F_rhoVyR, F_rhoVzR, F_BxR, F_ByR, F_BzR, F_EtotR] + = ideal_flux_vector_(rhoR, VxR, VyR, VzR, BxR, ByR, BzR, PR); + + // Non ideal contributions + hall_contribution_(rhoL, BxL, ByL, BzL, JxL, JyL, JzL, F_BxL, F_ByL, F_BzL, + F_EtotL); + + hall_contribution_(rhoR, BxR, ByR, BzR, JxR, JyR, JzR, F_BxR, F_ByR, F_BzR, + F_EtotR); + + resistive_contributions_(eta_, JxL, JyL, JzL, BxL, ByL, BzL, F_BxL, F_ByL, F_BzL, + F_EtotL); + + resistive_contributions_(eta_, JxR, JyR, JzR, BxR, ByR, BzR, F_BxR, F_ByR, F_BzR, + F_EtotR); + + auto [LaplJxL, LaplJxR] = reconstructed_laplacian(Jx, index); + auto [LaplJyL, LaplJyR] = reconstructed_laplacian(Jy, index); + auto [LaplJzL, LaplJzR] = reconstructed_laplacian(Jz, index); + + resistive_contributions_(nu_, LaplJxL, LaplJyL, LaplJzL, BxL, ByL, BzL, F_BxL, + F_ByL, F_BzL, F_EtotL); + + resistive_contributions_(nu_, LaplJxR, LaplJyR, LaplJzR, BxR, ByR, BzR, F_BxR, + F_ByR, F_BzR, F_EtotR); + + auto uL = std::forward_as_tuple(rhoL, VxL, VyL, VzL, BxL, ByL, BzL, PL); + auto uR = std::forward_as_tuple(rhoR, VxR, VyR, VzR, BxR, ByR, BzR, PR); + auto fL = std::forward_as_tuple(F_rhoL, F_rhoVxL, F_rhoVyL, F_rhoVzL, F_BxL, F_ByL, F_BzL, + F_EtotL); + auto fR = std::forward_as_tuple(F_rhoR, F_rhoVxR, F_rhoVyR, F_rhoVzR, F_BxR, F_ByR, F_BzR, + F_EtotR); + + auto [rho_, rhoVx_, rhoVy_, rhoVz_, Bx_, By_, Bz_, Etot_] + = riemann_solver_(uL, uR, fL, fR); + + if constexpr (dimension == 1) + { + rho_flux(index[0]) = rho_; + rhoV_flux(Component::X)(index[0]) = rhoVx_; + rhoV_flux(Component::Y)(index[0]) = rhoVy_; + rhoV_flux(Component::Z)(index[0]) = rhoVz_; + B_flux(Component::X)(index[0]) = Bx_; + B_flux(Component::Y)(index[0]) = By_; + B_flux(Component::Z)(index[0]) = Bz_; + Etot_flux(index[0]) = Etot_; + } + if constexpr (dimension == 2) + { + rho_flux(index[0], index[1]) = rho_; + rhoV_flux(Component::X)(index[0], index[1]) = rhoVx_; + rhoV_flux(Component::Y)(index[0], index[1]) = rhoVy_; + rhoV_flux(Component::Z)(index[0], index[1]) = rhoVz_; + B_flux(Component::X)(index[0], index[1]) = Bx_; + B_flux(Component::Y)(index[0], index[1]) = By_; + B_flux(Component::Z)(index[0], index[1]) = Bz_; + Etot_flux(index[0], index[1]) = Etot_; + } + if constexpr (dimension == 3) + { + rho_flux(index[0], index[1], index[2]) = rho_; + rhoV_flux(Component::X)(index[0], index[1], index[2]) = rhoVx_; + rhoV_flux(Component::Y)(index[0], index[1], index[2]) = rhoVy_; + rhoV_flux(Component::Z)(index[0], index[1], index[2]) = rhoVz_; + B_flux(Component::X)(index[0], index[1], index[2]) = Bx_; + B_flux(Component::Y)(index[0], index[1], index[2]) = By_; + B_flux(Component::Z)(index[0], index[1], index[2]) = Bz_; + Etot_flux(index[0], index[1], index[2]) = Etot_; + } + } + + template + auto riemann_solver_(auto const& uL, auto const& uR, auto const& fL, auto const& fR) const + { + auto const& [rhoL, VxL, VyL, VzL, BxL, ByL, BzL, PL] = uL; + auto const& [rhoR, VxR, VyR, VzR, BxR, ByR, BzR, PR] = uR; + auto const& [F_rhoL, F_rhoVxL, F_rhoVyL, F_rhoVzL, F_BxL, F_ByL, F_BzL, F_EtotL] = fL; + auto const& [F_rhoR, F_rhoVxR, F_rhoVyR, F_rhoVzR, F_BxR, F_ByR, F_BzR, F_EtotR] = fR; + + // Convert to conserved variables + auto rhoVxL = rhoL * VxL; + auto rhoVyL = rhoL * VyL; + auto rhoVzL = rhoL * VzL; + auto EtotL = PL / (gamma_ - 1) + 0.5 * rhoL * (VxL * VxL + VyL * VyL + VzL * VzL) + + 0.5 * (BxL * BxL + ByL * ByL + BzL * BzL); + + auto rhoVxR = rhoR * VxR; + auto rhoVyR = rhoR * VyR; + auto rhoVzR = rhoR * VzR; + auto EtotR = PR / (gamma_ - 1) + 0.5 * rhoR * (VxR * VxR + VyR * VyR + VzR * VzR) + + 0.5 * (BxR * BxR + ByR * ByR + BzR * BzR); + + auto uL_ = std::forward_as_tuple(rhoL, rhoVxL, rhoVyL, rhoVzL); + auto uR_ = std::forward_as_tuple(rhoR, rhoVxR, rhoVyR, rhoVzR); + auto fL_ = std::forward_as_tuple(F_rhoL, F_rhoVxL, F_rhoVyL, F_rhoVzL); + auto fR_ = std::forward_as_tuple(F_rhoR, F_rhoVxR, F_rhoVyR, F_rhoVzR); + + auto ubL = std::forward_as_tuple(BxL, ByL, BzL, EtotL); + auto ubR = std::forward_as_tuple(BxR, ByR, BzR, EtotR); + auto fbL = std::forward_as_tuple(F_BxL, F_ByL, F_BzL, F_EtotL); + auto fbR = std::forward_as_tuple(F_BxR, F_ByR, F_BzR, F_EtotR); + + // for rusanov riemann solver + auto const [S, Sb] = rusanov_speeds_(uL, uR); + auto [rho_, rhoVx_, rhoVy_, rhoVz_] = rusanov_(uL_, uR_, fL_, fR_, S); + auto [Bx_, By_, Bz_, Etot_] = rusanov_(ubL, ubR, fbL, fbR, Sb); + + return std::make_tuple(rho_, rhoVx_, rhoVy_, rhoVz_, Bx_, By_, Bz_, Etot_); + } + + template + auto rusanov_speeds_(auto const& uL, auto const& uR) const + { + auto const& [rhoL, VxL, VyL, VzL, BxL, ByL, BzL, PL] = uL; + auto const& [rhoR, VxR, VyR, VzR, BxR, ByR, BzR, PR] = uR; + auto BdotBL = BxL * BxL + ByL * ByL + BzL * BzL; + auto BdotBR = BxR * BxR + ByR * ByR + BzR * BzR; + + auto compute_speeds = [&](auto rhoL, auto rhoR, auto BdotBL, auto BdotBR, auto PL, auto PR, + auto Bcomp, auto Vcomp) { + auto cfastL = compute_fast_magnetosonic_(rhoL, Bcomp, BdotBL, PL); + auto cfastR = compute_fast_magnetosonic_(rhoR, Bcomp, BdotBR, PR); + auto cwL = compute_whistler_(layout_->inverseMeshSize(direction), rhoL, BdotBL); + auto cwR = compute_whistler_(layout_->inverseMeshSize(direction), rhoR, BdotBR); + auto S = std::max(std::abs(Vcomp) + cfastL, std::abs(Vcomp) + cfastR); + auto Sb = std::max(std::abs(Vcomp) + cfastL + cwL, std::abs(Vcomp) + cfastR + cwR); + return std::make_tuple(S, Sb); + }; + + if constexpr (direction == Direction::X) + return compute_speeds(rhoL, rhoR, BdotBL, BdotBR, PL, PR, BxL, VxL); + else if constexpr (direction == Direction::Y) + return compute_speeds(rhoL, rhoR, BdotBL, BdotBR, PL, PR, ByL, VyL); + else if constexpr (direction == Direction::Z) + return compute_speeds(rhoL, rhoR, BdotBL, BdotBR, PL, PR, BzL, VzL); + } + + auto rusanov_(auto const& uL, auto const& uR, auto const& fL, auto const& fR, + auto const S) const + { + // to be used 2 times in hall mhd (the second time for B with whisler contribution). + + auto constexpr N_elements = std::tuple_size_v>; + + return for_N([&](auto i) { + return (std::get(fL) + std::get(fR)) * 0.5 + - S * (std::get(uR) - std::get(uL)) * 0.5; + }); + } + + auto compute_fast_magnetosonic_(auto const& rho, auto const& B, auto const& BdotB, + auto const& P) const + { + auto Sound = std::sqrt((gamma_ * P) / rho); + auto AlfvenDir = std::sqrt(B * B / rho); // diectionnal alfven + auto Alfven = std::sqrt(BdotB / rho); + + auto c02 = Sound * Sound; + auto cA2 = Alfven * Alfven; + auto cAdir2 = AlfvenDir * AlfvenDir; + + return std::sqrt((c02 + cA2) * 0.5 + + std::sqrt((c02 + cA2) * (c02 + cA2) - 4.0 * c02 * cAdir2)); + } + + auto compute_whistler_(auto const& invMeshSize, auto const& rho, auto const& BdotB) const + { + auto vw = std::sqrt(1 + 0.25 * invMeshSize * invMeshSize) + 0.5 * invMeshSize; + return std::sqrt(BdotB) * vw / rho; + } + + template + auto ideal_flux_vector_(auto const& rho, auto const& Vx, auto const& Vy, auto const& Vz, + auto const& Bx, auto const& By, auto const& Bz, auto const& P) const + { + auto GeneralisedPressure = P + 0.5 * (Bx * Bx + By * By + Bz * Bz); + auto TotalEnergy = P / (gamma_ - 1) + 0.5 * rho * (Vx * Vx + Vy * Vy + Vz * Vz) + + 0.5 * (Bx * Bx + By * By + Bz * Bz); + if constexpr (direction == Direction::X) + { + auto F_rho = rho * Vx; + auto F_rhoVx = rho * Vx * Vx + GeneralisedPressure + Bx * Bx; + auto F_rhoVy = rho * Vx * Vy + Bx * By; + auto F_rhoVz = rho * Vx * Vz + Bx * Bz; + auto F_Bx = 0.0; + auto F_By = By * Vx - Vy * Bx; + auto F_Bz = Bz * Vx - Vz * Bx; + auto F_Etot + = (TotalEnergy + GeneralisedPressure) * Vx - Bx * (Vx * Bx + Vy * By + Vz * Bz); + + return std::make_tuple(F_rho, F_rhoVx, F_rhoVy, F_rhoVz, F_Bx, F_By, F_Bz, F_Etot); + } + if constexpr (direction == Direction::Y) + { + auto F_rho = rho * Vy; + auto F_rhoVx = rho * Vy * Vx + By * Bx; + auto F_rhoVy = rho * Vy * Vy + GeneralisedPressure + By * By; + auto F_rhoVz = rho * Vy * Vz + By * Bz; + auto F_Bx = Bx * Vy - Vx * By; + auto F_By = 0.0; + auto F_Bz = Bz * Vy - Vz * By; + auto F_Etot + = (TotalEnergy + GeneralisedPressure) * Vy - By * (Vx * Bx + Vy * By + Vz * Bz); + + return std::make_tuple(F_rho, F_rhoVx, F_rhoVy, F_rhoVz, F_Bx, F_By, F_Bz, F_Etot); + } + if constexpr (direction == Direction::Z) + { + auto F_rho = rho * Vz; + auto F_rhoVx = rho * Vz * Vx + Bz * Bx; + auto F_rhoVy = rho * Vz * Vy + Bz * By; + auto F_rhoVz = rho * Vz * Vz + GeneralisedPressure + Bz * Bz; + auto F_Bx = Bx * Vz - Vx * Bz; + auto F_By = By * Vz - Vy * Bz; + auto F_Bz = 0.0; + auto F_Etot + = (TotalEnergy + GeneralisedPressure) * Vz - Bz * (Vx * Bx + Vy * By + Vz * Bz); + + return std::make_tuple(F_rho, F_rhoVx, F_rhoVy, F_rhoVz, F_Bx, F_By, F_Bz, F_Etot); + } + } + + template + void hall_contribution_(auto const& rho, auto const& Bx, auto const& By, auto const& Bz, + auto const& Jx, auto const& Jy, auto const& Jz, auto& F_Bx, auto& F_By, + auto& F_Bz, auto& F_Etot) const + { + auto invRho = 1.0 / rho; + + auto JxB_x = Jy * Bz - Jz * By; + auto JxB_y = Jz * Bx - Jx * Bz; + auto JxB_z = Jx * By - Jy * Bx; + + auto BdotJ = Bx * Jx + By * Jy + Bz * Jz; + auto BdotB = Bx * Bx + By * By + Bz * Bz; + + if constexpr (direction == Direction::X) + { + F_By += -JxB_z * invRho; + F_Bz += JxB_y * invRho; + F_Etot += (BdotJ * Bx - BdotB * Jx) * invRho; + } + if constexpr (direction == Direction::Y) + { + F_Bx += JxB_z * invRho; + F_Bz += -JxB_x * invRho; + F_Etot += (BdotJ * By - BdotB * Jy) * invRho; + } + if constexpr (direction == Direction::Z) + { + F_Bx += -JxB_y * invRho; + F_By += JxB_x * invRho; + F_Etot += (BdotJ * Bz - BdotB * Jz) * invRho; + } + } + + template + void resistive_contributions_(auto const& pc, auto const& Jx, auto const& Jy, auto const& Jz, + auto const& Bx, auto const& By, auto const& Bz, auto& F_Bx, + auto& F_By, auto& F_Bz, auto& F_Etot) const + // Can be used for both resistivity with J and eta and hyper resistivity with laplJ and nu + { + if constexpr (direction == Direction::X) + { + F_By += -Jz * pc; + F_Bz += Jy * pc; + F_Etot += (Jy * Bz - Jz * By) * pc; + } + if constexpr (direction == Direction::Y) + { + F_Bx += Jz * pc; + F_Bz += -Jx * pc; + F_Etot += (Jz * Bx - Jx * Bz) * pc; + } + if constexpr (direction == Direction::Y) + { + F_Bx += -Jy * pc; + F_By += Jx * pc; + F_Etot += (Jx * By - Jy * Bx) * pc; + } + } + + template + auto reconstructed_laplacian(Field const& J, MeshIndex index) const + { + auto d2 + = [&](Direction dir, auto const& prevValue, auto const& Value, auto const& nextValue) { + return (layout_->inverseMeshSize(dir)) * (layout_->inverseMeshSize(dir)) + * (prevValue - 2.0 * Value + nextValue); + }; + + auto JL = reconstruct_uL_(J, index); + auto JR = reconstruct_uR_(J, index); + + if constexpr (dimension == 1) + { + MeshIndex<1> prevX = index; + prevX[0] -= 1; + MeshIndex<1> nextX = index; + nextX[0] += 1; + + auto JL_X_1 = reconstruct_uL_(J, prevX); + auto JR_X_1 = reconstruct_uR_(J, prevX); + auto JL_X1 = reconstruct_uL_(J, nextX); + auto JR_X1 = reconstruct_uR_(J, nextX); + + auto LaplJL = d2(Direction::X, JL_X_1, JL, JL_X1); + auto LaplJR = d2(Direction::X, JR_X_1, JR, JR_X1); + + return std::make_tuple(LaplJL, LaplJR); + } + if constexpr (dimension == 2) + { + MeshIndex<2> prevX = index; + prevX[0] -= 1; + MeshIndex<2> nextX = index; + nextX[0] += 1; + + MeshIndex<2> prevY = index; + prevY[1] -= 1; + MeshIndex<2> nextY = index; + nextY[1] += 1; + + auto JL_X_1 = reconstruct_uL_(J, prevX); + auto JR_X_1 = reconstruct_uR_(J, prevX); + auto JL_X1 = reconstruct_uL_(J, nextX); + auto JR_X1 = reconstruct_uR_(J, nextX); + + auto JL_Y_1 = reconstruct_uL_(J, prevY); + auto JR_Y_1 = reconstruct_uR_(J, prevY); + auto JL_Y1 = reconstruct_uL_(J, nextY); + auto JR_Y1 = reconstruct_uR_(J, nextY); + + auto LaplJL = d2(Direction::X, JL_X_1, JL, JL_X1) + d2(Direction::Y, JL_Y_1, JL, JL_Y1); + auto LaplJR = d2(Direction::X, JR_X_1, JR, JR_X1) + d2(Direction::Y, JR_Y_1, JR, JR_Y1); + + return std::make_tuple(LaplJL, LaplJR); + } + if constexpr (dimension == 3) + { + MeshIndex<3> prevX = index; + prevX[0] -= 1; + MeshIndex<3> nextX = index; + nextX[0] += 1; + + MeshIndex<3> prevY = index; + prevY[1] -= 1; + MeshIndex<3> nextY = index; + nextY[1] += 1; + + MeshIndex<3> prevZ = index; + prevZ[2] -= 1; + MeshIndex<3> nextZ = index; + nextZ[2] += 1; + + auto JL_X_1 = reconstruct_uL_(J, prevX); + auto JR_X_1 = reconstruct_uR_(J, prevX); + auto JL_X1 = reconstruct_uL_(J, nextX); + auto JR_X1 = reconstruct_uR_(J, nextX); + + auto JL_Y_1 = reconstruct_uL_(J, prevY); + auto JR_Y_1 = reconstruct_uR_(J, prevY); + auto JL_Y1 = reconstruct_uL_(J, nextY); + auto JR_Y1 = reconstruct_uR_(J, nextY); + + auto JL_Z_1 = reconstruct_uL_(J, prevZ); + auto JR_Z_1 = reconstruct_uR_(J, prevZ); + auto JL_Z1 = reconstruct_uL_(J, nextZ); + auto JR_Z1 = reconstruct_uR_(J, nextZ); + + auto LaplJL = d2(Direction::X, JL_X_1, JL, JL_X1) + d2(Direction::Y, JL_Y_1, JL, JL_Y1) + + d2(Direction::Z, JL_Z_1, JL, JL_Z1); + auto LaplJR = d2(Direction::X, JR_X_1, JR, JR_X1) + d2(Direction::Y, JR_Y_1, JR, JR_Y1) + + d2(Direction::Z, JR_Z_1, JR, JR_Z1); + + return std::make_tuple(LaplJL, LaplJR); + } + } + + template + auto reconstruct_uL_(Field const& F, MeshIndex index) const + { + return constant_uL_(F, index); + } + + template + auto reconstruct_uR_(Field const& F, MeshIndex index) const + { + return constant_uR_(F, index); + } + + template + auto constant_uL_(Field const& F, MeshIndex index) const + { + auto fieldCentering = layout_->centering(F.physicalQuantity()); + + if constexpr (dimension == 1) + { + return F(layout_->prevIndex(fieldCentering[dirX], index[0])); + } + else if constexpr (dimension == 2) + { + if constexpr (direction == Direction::X) + { + return F(layout_->prevIndex(fieldCentering[dirX], index[0]), index[1]); + } + else if constexpr (direction == Direction::Y) + { + return F(index[0], layout_->prevIndex(fieldCentering[dirY], index[1])); + } + } + else if constexpr (dimension == 3) + { + if constexpr (direction == Direction::X) + { + return F(layout_->prevIndex(fieldCentering[dirX], index[0]), index[1], index[2]); + } + else if constexpr (direction == Direction::Y) + { + return F(index[0], layout_->prevIndex(fieldCentering[dirY], index[1]), index[2]); + } + else if constexpr (direction == Direction::Z) + { + return F(index[0], index[1], layout_->prevIndex(fieldCentering[dirY], index[2])); + } + } + } + + template + auto constant_uR_(Field const& F, MeshIndex index) const + { + if constexpr (dimension == 1) + { + return F(index[0]); + } + else if constexpr (dimension == 2) + { + return F(index[0], index[1]); + } + else if constexpr (dimension == 3) + { + return F(index[0], index[1], index[2]); + } + } +}; + +} // namespace PHARE::core + +#endif diff --git a/src/core/numerics/mhd_reconstruction/mhd_reconstruction.hpp b/src/core/numerics/mhd_reconstruction/mhd_reconstruction.hpp deleted file mode 100644 index bf35dde7d..000000000 --- a/src/core/numerics/mhd_reconstruction/mhd_reconstruction.hpp +++ /dev/null @@ -1,84 +0,0 @@ -#ifndef PHARE_CORE_NUMERICS_MHD_RECONSTRUCTION_HPP -#define PHARE_CORE_NUMERICS_MHD_RECONSTRUCTION_HPP - -#include "core/data/grid/gridlayout.hpp" -#include "core/data/grid/gridlayout_utils.hpp" -#include "core/data/grid/gridlayoutdefs.hpp" -#include "core/utilities/index/index.hpp" -#include - -namespace PHARE::core -{ -template -class Reconstruction : public LayoutHolder -{ - constexpr static auto dimension = GridLayout::dimension; - using LayoutHolder::layout_; - -public: - template - void operator()(Field const& F, Field& uL, Field& uR) - { - if (!this->hasLayout()) - throw std::runtime_error("Error - Reconstruction - GridLayout not set, cannot proceed " - "to reconstruction"); - - if (!(F.isUsable() && uL.isUsable() && uR.isUsable())) - throw std::runtime_error( - "Error - Reconstruction - not all Field parameters are usable"); - - - layout_->evalOnBox( - uL, [&](auto&... args) { this->template reconstruct_uL_(F, uL, args...); }); - - layout_->evalOnBox( - uR, [&](auto&... args) { this->template reconstruct_uL_(F, uR, args...); }); - } - -private: - template - void reconstruct_uL_(Field const& F, Field& uL, Indexes const&... ijk) const - { - uL(ijk...) = this->template constant_uL_(F, {ijk...}); - } - - template - void reconstruct_uR_(Field const& F, Field& uR, Indexes const&... ijk) const - { - uR(ijk...) = this->template constant_uR_(F, {ijk...}); - } - - template - auto constant_uL_(Field const& F, MeshIndex index) const - { - auto fieldCentering = layout_->centering(F.physicalQuantity()); - - std::size_t dir; - if (direction == Direction::X) - dir = PHARE::core::dirX; - else if (direction == Direction::Y) - dir = PHARE::core::dirY; - else - dir = PHARE::core::dirZ; - - return (layout_->prevIndex(fieldCentering[dir], index[0])); - } - - template - auto constant_uR_(Field const& F, MeshIndex index) const - { - std::size_t dir; - if (direction == Direction::X) - dir = PHARE::core::dirX; - else if (direction == Direction::Y) - dir = PHARE::core::dirY; - else - dir = PHARE::core::dirZ; - - return F(index[0]); - } -}; - -} // namespace PHARE::core - -#endif diff --git a/tests/core/data/field/test_field_fixtures_mhd.hpp b/tests/core/data/field/test_field_fixtures_mhd.hpp index 69f4c5044..4b4a98f37 100644 --- a/tests/core/data/field/test_field_fixtures_mhd.hpp +++ b/tests/core/data/field/test_field_fixtures_mhd.hpp @@ -1,5 +1,5 @@ -#ifndef PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP -#define PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_HPP +#ifndef PHARE_TEST_CORE_DATA_TEST_FIELD_MHD_HPP +#define PHARE_TEST_CORE_DATA_TEST_FIELD_MHD_HPP #include "core/data/field/field.hpp" #include "core/mhd/mhd_quantities.hpp" diff --git a/tests/core/data/field/test_usable_field_fixtures_mhd.hpp b/tests/core/data/field/test_usable_field_fixtures_mhd.hpp new file mode 100644 index 000000000..86da4b580 --- /dev/null +++ b/tests/core/data/field/test_usable_field_fixtures_mhd.hpp @@ -0,0 +1,39 @@ +#ifndef PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_MHD_HPP +#define PHARE_TEST_CORE_DATA_TEST_FIELD_FIXTURES_MHD_HPP + + +#include "core/mhd/mhd_quantities.hpp" +#include "core/data/field/field.hpp" +#include "core/data/grid/grid.hpp" + +namespace PHARE::core +{ + +template +class UsableFieldMHD : Field +{ +public: + auto static constexpr dimension = dim; + using Super = Field; + using Grid_t = Grid, MHDQuantity::Scalar>; + using scalar_t = MHDQuantity::Scalar; + + template + UsableFieldMHD(std::string const& name, GridLayout const& layout, scalar_t qty) + : Super{name, qty} + , xyz{name, layout, qty} + { + super().setBuffer(&xyz); + } + void set_on(Super& field) { field.setBuffer(&xyz); } + + Super& super() { return *this; } + Super const& super() const { return *this; } + +protected: + Grid_t xyz; +}; + +} // namespace PHARE::core + +#endif diff --git a/tests/core/data/tensorfield/test_tensorfield_fixtures_mhd.hpp b/tests/core/data/tensorfield/test_tensorfield_fixtures_mhd.hpp index f49d85900..b75699a9c 100644 --- a/tests/core/data/tensorfield/test_tensorfield_fixtures_mhd.hpp +++ b/tests/core/data/tensorfield/test_tensorfield_fixtures_mhd.hpp @@ -12,8 +12,8 @@ namespace PHARE::core { /* -A UsableTensorFieldMHD is an extension of the TensorField view that owns memory for components and sets -the view pointers. It is useful for tests to easily declare usable (== set views) tensors +A UsableTensorFieldMHD is an extension of the TensorField view that owns memory for components and +sets the view pointers. It is useful for tests to easily declare usable (== set views) tensors */ template class UsableTensorFieldMHD : public TensorField, MHDQuantity, rank_> @@ -43,16 +43,15 @@ class UsableTensorFieldMHD : public TensorField, MHDQuantity, rank } Super& super() { return *this; } - Super& super() const { return *this; } + Super const& super() const { return *this; } protected: template auto static make_grids(ComponentNames const& compNames, GridLayout const& layout, tensor_t qty) { auto qts = MHDQuantity::componentsQuantities(qty); - return for_N([&](auto i) { - return Grid_t{compNames[i], qts[i], layout.allocSize(qts[i])}; - }); + return for_N( + [&](auto i) { return Grid_t{compNames[i], qts[i], layout.allocSize(qts[i])}; }); } std::array xyz; diff --git a/tests/core/numerics/mhd_solver/test_mhd_solver.cpp b/tests/core/numerics/mhd_solver/test_mhd_solver.cpp index 0a5a288a9..1fbf1c00f 100644 --- a/tests/core/numerics/mhd_solver/test_mhd_solver.cpp +++ b/tests/core/numerics/mhd_solver/test_mhd_solver.cpp @@ -1,4 +1,6 @@ #include +#include +#include #include #include #include @@ -12,11 +14,15 @@ #include "core/data/grid/gridlayoutimplyee_mhd.hpp" #include "core/numerics/boundary_condition/boundary_condition.hpp" #include "phare_core.hpp" +#include "tests/core/data/field/test_field_fixtures_mhd.hpp" #include "tests/core/data/gridlayout/test_gridlayout.hpp" #include "tests/core/data/mhd_state/test_mhd_state_fixtures.hpp" #include "amr/solvers/solver_mhd.hpp" #include "amr/solvers/solver.hpp" +#include "tests/core/data/field/test_usable_field_fixtures_mhd.hpp" +#include "tests/core/data/vecfield/test_vecfield_fixtures_mhd.hpp" + constexpr std::uint32_t cells = 65; constexpr std::size_t dim = 1, interp = 1; @@ -24,23 +30,130 @@ using YeeLayout_t = PHARE::core::GridLayoutImplYeeMHD; using GridLayoutMHD = PHARE::core::GridLayout; using GridLayout_t = TestGridLayout; using BoundaryCondition = PHARE::core::BoundaryCondition; +using MHDQuantity = PHARE::core::MHDQuantity; +using FieldMHD = PHARE::core::FieldMHD; +using VecFieldMHD = PHARE::core::VecField; + +using namespace PHARE::initializer; -using DummyModelView = PHARE::core::UsableMHDState; +PHAREDict getDict() +{ + PHAREDict dict; + dict["godunov"]["resistivity"] = 0.01; + dict["godunov"]["hyper_resistivity"] = 0.01; + dict["godunov"]["heat_capacity_ratio"] = 0.01; + + return dict; +} -template -class UsableMHDStateWrapper : public PHARE::solver::ISolverModelView +struct DummyModelViewConstructor : public PHARE::solver::ISolverModelView { -public: - UsableMHDStateWrapper(PHARE::core::UsableMHDState& mhdState) - : mhdState_(mhdState) + using GodunovFluxes_t = PHARE::core::GodunovFluxes; + + DummyModelViewConstructor(GridLayout_t const& layout) + : rho{"rho", layout, MHDQuantity::Scalar::rho} + , V{"V", layout, MHDQuantity::Vector::V} + , B_CT{"B_CT", layout, MHDQuantity::Vector::B_CT} + , P{"P", layout, MHDQuantity::Scalar::P} + , J{"J", layout, MHDQuantity::Vector::J} + + , rho_x{"rho_x", layout, MHDQuantity::Scalar::ScalarFlux_x} + , rhoV_x{"V_x", layout, MHDQuantity::Vector::VecFlux_x} + , B_x{"B_x", layout, MHDQuantity::Vector::VecFlux_x} + , Etot_x{"Etot_x", layout, MHDQuantity::Scalar::ScalarFlux_x} + + , rho_y{"rho_y", layout, MHDQuantity::Scalar::ScalarFlux_y} + , rhoV_y{"V_y", layout, MHDQuantity::Vector::VecFlux_y} + , B_y{"B_y", layout, MHDQuantity::Vector::VecFlux_y} + , Etot_y{"Etot_y", layout, MHDQuantity::Scalar::ScalarFlux_y} + + , rho_z{"rho_z", layout, MHDQuantity::Scalar::ScalarFlux_z} + , rhoV_z{"V_z", layout, MHDQuantity::Vector::VecFlux_z} + , B_z{"B_z", layout, MHDQuantity::Vector::VecFlux_z} + , Etot_z{"Etot_z", layout, MHDQuantity::Scalar::ScalarFlux_z} + + , layouts{layout} { } -private: - PHARE::core::UsableMHDState& mhdState_; // Reference to the original state + PHARE::core::UsableFieldMHD rho; + PHARE::core::UsableVecFieldMHD V; + PHARE::core::UsableVecFieldMHD B_CT; + PHARE::core::UsableFieldMHD P; + PHARE::core::UsableVecFieldMHD J; + + PHARE::core::UsableFieldMHD rho_x; + PHARE::core::UsableVecFieldMHD rhoV_x; + PHARE::core::UsableVecFieldMHD B_x; + PHARE::core::UsableFieldMHD Etot_x; + + PHARE::core::UsableFieldMHD rho_y; + PHARE::core::UsableVecFieldMHD rhoV_y; + PHARE::core::UsableVecFieldMHD B_y; + PHARE::core::UsableFieldMHD Etot_y; + + PHARE::core::UsableFieldMHD rho_z; + PHARE::core::UsableVecFieldMHD rhoV_z; + PHARE::core::UsableVecFieldMHD B_z; + PHARE::core::UsableFieldMHD Etot_z; + + GridLayout_t layouts; }; +struct DummyModelView : public PHARE::solver::ISolverModelView +{ + using GodunovFluxes_t = PHARE::core::GodunovFluxes; + + DummyModelView(DummyModelViewConstructor& construct) + { + rho.push_back(&construct.rho.super()); + V.push_back(&construct.V.super()); + B_CT.push_back(&construct.B_CT.super()); + P.push_back(&construct.P.super()); + J.push_back(&construct.J.super()); + + rho_x.push_back(&construct.rho_x.super()); + rhoV_x.push_back(&construct.rhoV_x.super()); + B_x.push_back(&construct.B_x.super()); + Etot_x.push_back(&construct.Etot_x.super()); + + rho_y.push_back(&construct.rho_y.super()); + rhoV_y.push_back(&construct.rhoV_y.super()); + B_y.push_back(&construct.B_y.super()); + Etot_y.push_back(&construct.Etot_y.super()); + + rho_z.push_back(&construct.rho_z.super()); + rhoV_z.push_back(&construct.rhoV_z.super()); + B_z.push_back(&construct.B_z.super()); + Etot_z.push_back(&construct.Etot_z.super()); + + layouts.push_back(&construct.layouts); + } + + std::vector rho; + std::vector V; + std::vector B_CT; + std::vector P; + std::vector J; + + std::vector rho_x; + std::vector rhoV_x; + std::vector B_x; + std::vector Etot_x; + + std::vector rho_y; + std::vector rhoV_y; + std::vector B_y; + std::vector Etot_y; + + std::vector rho_z; + std::vector rhoV_z; + std::vector B_z; + std::vector Etot_z; + + std::vector layouts; +}; struct DummyHierarchy { @@ -68,12 +181,11 @@ struct DummyRessourcesManager struct DummyMHDModel : public PHARE::solver::IPhysicalModel { - using field_type - = PHARE::core::FieldMHD; // typename PHARE::core::UsableMHDState::Grid_t; + using field_type = FieldMHD; + using vecfield_type = VecFieldMHD; using gridlayout_type = GridLayout_t; static constexpr std::size_t dimension = dim; static constexpr auto model_name = "mhd_model"; - DummyRessourcesManager ressourcesManager; }; @@ -150,11 +262,11 @@ class DummyMessenger : public PHARE::amr::IMessenger - TestMHDSolver; + TestMHDSolver(getDict()); - UsableMHDStateWrapper dummy_view(state); // torm DummyMessenger dummy_messenger;