Skip to content

Commit

Permalink
initial branch commit: new architecture for reconstruction + flux com…
Browse files Browse the repository at this point in the history 
…putation

improved test, fixed some bugs in godunov-fluxes

better timing directory creation (PHAREHUB#914)

safer nu (PHAREHUB#911)

try fallback for dl on keyerror + ruffage (PHAREHUB#916)

convenient utils for mpi/hierarchies

keep pyattrs on compute_from_hier

nan min/max to handle possible nan ghosts (PHAREHUB#923)

* nan min/max to handle possible nan ghosts

* flatten

rm atefact file

fixed missing template keyword for macos-12 build

more explicit unwrapping in godunov fluxes with variadic arguments of unclear size

fixed lambda capture problem

fixed lambda capture problem
  • Loading branch information
@UCaromel
UCaromel committed Jan 10, 2025
1 parent b4528b1 commit e1574ef
Show file tree
Hide file tree
Showing 13 changed files with 1,133 additions and 363 deletions.
127 changes: 52 additions & 75 deletions src/amr/solvers/solver_mhd.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand All @@ -31,27 +32,41 @@ class SolverMHD : public ISolver<AMR_Types>
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<AMR_Types>;
using IMessenger = amr::IMessenger<IPhysicalModel_t>;
using Direction = core::Direction;

using Reconstruction_t = core::Reconstruction<GridLayout>;
using GodunovFluxes_t = typename MHDModelView<MHDModel>::GodunovFluxes_t;
using Ampere_t = typename MHDModelView<MHDModel>::Ampere_t;

std::vector<field_type> uL_x;
std::vector<field_type> uR_x;
std::vector<field_type> uL_y;
std::vector<field_type> uR_y;
std::vector<field_type> uL_z;
std::vector<field_type> 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<AMR_Types>{"MHDSolver"}
, godunov_{dict["godunov"]}
{
}

Expand Down Expand Up @@ -82,10 +97,7 @@ class SolverMHD : public ISolver<AMR_Types>
}

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,
Expand Down Expand Up @@ -123,9 +135,7 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::advanceLevel(
auto& fromCoarser = dynamic_cast<Messenger&>(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);

Expand All @@ -135,68 +145,34 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::advanceLevel(
}

template<typename MHDModel, typename AMR_Types, typename Messenger, typename ModelViews_t>
void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::reconstruction_(
void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::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::size_t... I>(std::index_sequence<I...>) {
((reconstruct_.template operator()<Direction::X>(std::get<I>(Fields), uL_x[I],
uR_x[I])),
...);
}(std::make_index_sequence<8>{});
}
if constexpr (dimension == 2)
{
[&]<std::size_t... I>(std::index_sequence<I...>) {
((reconstruct_.template operator()<Direction::X>(std::get<I>(Fields), uL_x[I],
uR_x[I])),
...);
((reconstruct_.template operator()<Direction::Y>(std::get<I>(Fields), uL_y[I],
uR_y[I])),
...);
}(std::make_index_sequence<8>{});
}
if constexpr (dimension == 3)
{
[&]<std::size_t... I>(std::index_sequence<I...>) {
((reconstruct_.template operator()<Direction::X>(std::get<I>(Fields), uL_x[I],
uR_x[I])),
...);
((reconstruct_.template operator()<Direction::Y>(std::get<I>(Fields), uL_y[I],
uR_y[I])),
...);
((reconstruct_.template operator()<Direction::Z>(std::get<I>(Fields), uL_z[I],
uR_z[I])),
...);
}(std::make_index_sequence<8>{});
}
},
Fields);
}
PHARE_LOG_SCOPE(1, "SolverMHD::godunov_fluxes_");

template<typename MHDModel, typename AMR_Types, typename Messenger, typename ModelViews_t>
void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::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<typename MHDModel, typename AMR_Types, typename Messenger, typename ModelViews_t>
void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::FV_cycle_(level_t& level,
ModelViews_t& views,
Expand All @@ -221,6 +197,7 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::constrained_transp

auto dt = newTime - currentTime;

// Fill ghost for B_RS
// averaging B_RS(x, y, z)
// faraday
}
Expand Down
41 changes: 36 additions & 5 deletions src/amr/solvers/solver_mhd_model_view.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<typename GridLayout>
class GodunovFluxesTransformer
{
using core_type = PHARE::core::GodunovFluxes<GridLayout>;

template <typename MHDModel_>
class MHDModelView : public ISolverModelView {};
public:
template<typename Layout, typename Field, typename VecField, typename... Fluxes>
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<typename MHDModel_>
class MHDModelView : public ISolverModelView
{
public:
using MHDModel_t = MHDModel_;
using GridLayout = typename MHDModel_t::gridlayout_type;
using GodunovFluxes_t = GodunovFluxesTransformer<GridLayout>;
using Ampere_t = AmpereTransformer<GridLayout>;
};

}; // namespace PHARE::solver

#endif // PHARE_SOLVER_SOLVER_MHD_MODEL_VIEW_HPP
24 changes: 13 additions & 11 deletions src/core/data/grid/gridlayout.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::uint32_t>(quantity);
Expand Down Expand Up @@ -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<std::uint32_t>(quantity);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<std::uint32_t>(quantity);
Expand Down Expand Up @@ -809,7 +810,7 @@ namespace core
* @brief returns the centering of a scalar hybrid quantity in each directions
*/
NO_DISCARD constexpr static std::array<QtyCentering, dimension>
centering(Quantity::Scalar quantity)
centering(typename Quantity::Scalar quantity)
{
return GridLayoutImpl::centering(quantity);
}
Expand All @@ -818,7 +819,7 @@ namespace core
* @brief returns the centering of a vector hybrid quantity in each directions
*/
NO_DISCARD constexpr static std::array<std::array<QtyCentering, dimension>, 3>
centering(Quantity::Vector quantity)
centering(typename Quantity::Vector quantity)
{
return GridLayoutImpl::centering(quantity);
}
Expand All @@ -828,7 +829,8 @@ namespace core
* @return An std::array<std::uint32_t, dim> object, containing the size to which allocate
* arrays of an Quantity::Quantity 'qty' in every directions.
*/
NO_DISCARD std::array<std::uint32_t, dimension> allocSize(Quantity::Scalar qty) const
NO_DISCARD std::array<std::uint32_t, dimension>
allocSize(typename Quantity::Scalar qty) const
{
std::uint32_t iQty = static_cast<std::uint32_t>(qty);

Expand All @@ -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<std::uint32_t, dimension> allocSizeDerived(Quantity::Scalar qty,
Direction dir) const
NO_DISCARD std::array<std::uint32_t, dimension>
allocSizeDerived(typename Quantity::Scalar qty, Direction dir) const
{
std::uint32_t iDerivedDir = static_cast<std::uint32_t>(dir);
std::uint32_t iQty = static_cast<std::uint32_t>(qty);
Expand Down Expand Up @@ -893,7 +895,7 @@ namespace core
* ghost nodes
*/
NO_DISCARD std::array<std::uint32_t, dimension>
nbrPhysicalNodes(Quantity::Scalar hybQty) const
nbrPhysicalNodes(typename Quantity::Scalar hybQty) const
{
std::array<QtyCentering, dimension> centerings;

Expand All @@ -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<std::uint32_t>(qty);
std::uint32_t idir = static_cast<std::uint32_t>(dir);
Expand Down
47 changes: 27 additions & 20 deletions src/core/data/grid/gridlayoutdefs.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,26 +141,33 @@ namespace core
static constexpr std::uint32_t iJy = static_cast<std::uint32_t>(MHDQuantity::Scalar::Jy);
static constexpr std::uint32_t iJz = static_cast<std::uint32_t>(MHDQuantity::Scalar::Jz);

static constexpr std::uint32_t iBx_RSx
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bx_RSx);
static constexpr std::uint32_t iBy_RSx
= static_cast<std::uint32_t>(MHDQuantity::Scalar::By_RSx);
static constexpr std::uint32_t iBz_RSx
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bz_RSx);

static constexpr std::uint32_t iBx_RSy
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bx_RSy);
static constexpr std::uint32_t iBy_RSy
= static_cast<std::uint32_t>(MHDQuantity::Scalar::By_RSy);
static constexpr std::uint32_t iBz_RSy
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bz_RSy);

static constexpr std::uint32_t iBx_RSz
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bx_RSz);
static constexpr std::uint32_t iBy_RSz
= static_cast<std::uint32_t>(MHDQuantity::Scalar::By_RSz);
static constexpr std::uint32_t iBz_RSz
= static_cast<std::uint32_t>(MHDQuantity::Scalar::Bz_RSz);
static constexpr std::uint32_t iScalarFlux_x
= static_cast<std::uint32_t>(MHDQuantity::Scalar::ScalarFlux_x);
static constexpr std::uint32_t iScalarFlux_y
= static_cast<std::uint32_t>(MHDQuantity::Scalar::ScalarFlux_y);
static constexpr std::uint32_t iScalarFlux_z
= static_cast<std::uint32_t>(MHDQuantity::Scalar::ScalarFlux_z);

static constexpr std::uint32_t iVecFluxX_x
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxX_x);
static constexpr std::uint32_t iVecFluxY_x
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxY_x);
static constexpr std::uint32_t iVecFluxZ_x
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxZ_x);

static constexpr std::uint32_t iVecFluxX_y
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxX_y);
static constexpr std::uint32_t iVecFluxY_y
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxY_y);
static constexpr std::uint32_t iVecFluxZ_y
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxZ_y);

static constexpr std::uint32_t iVecFluxX_z
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxX_z);
static constexpr std::uint32_t iVecFluxY_z
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxY_z);
static constexpr std::uint32_t iVecFluxZ_z
= static_cast<std::uint32_t>(MHDQuantity::Scalar::VecFluxZ_z);
};
} // namespace core
} // namespace PHARE
Expand Down
Loading

0 comments on commit e1574ef

Please sign in to comment.