time integration refactor, added tvdrk2 and tvdrk3

Author: UCaromel

src/amr/solvers/solver_mhd.hpp

@@ -47,9 +47,7 @@ class SolverMHD : public ISolver<AMR_Types>
     using GodunovFluxes_t           = typename MHDModelView<MHDModel>::GodunovFluxes_t;
     using ToPrimitiveConverter_t    = typename MHDModelView<MHDModel>::ToPrimitiveConverter_t;
     using ToConservativeConverter_t = typename MHDModelView<MHDModel>::ToConservativeConverter_t;
-    using FiniteVolumeEuler_t       = typename MHDModelView<MHDModel>::FiniteVolumeEuler_t;
-    using ConstrainedTransport_t    = typename MHDModelView<MHDModel>::ConstrainedTransport_t;
-    using Faraday_t                 = typename MHDModelView<MHDModel>::Faraday_t;
+    using TimeIntegrator_t          = typename MHDModelView<MHDModel>::TimeIntegrator_t;
 
 
     // Flux calculations
@@ -72,20 +70,18 @@ class SolverMHD : public ISolver<AMR_Types>
     FieldT rho1{"rho1", MHDQuantity::Scalar::rho};
     VecFieldT rhoV1{"rhoV1", MHDQuantity::Vector::rhoV};
     VecFieldT B1{"B1", MHDQuantity::Vector::B};
-    FieldT Etot1{"rho1", MHDQuantity::Scalar::Etot};
+    FieldT Etot1{"Etot1", MHDQuantity::Scalar::Etot};
 
     FieldT rho2{"rho2", MHDQuantity::Scalar::rho};
     VecFieldT rhoV2{"rhoV2", MHDQuantity::Vector::rhoV};
     VecFieldT B2{"B2", MHDQuantity::Vector::B};
-    FieldT Etot2{"rho2", MHDQuantity::Scalar::Etot};
+    FieldT Etot2{"Etot2", MHDQuantity::Scalar::Etot};
 
     Ampere_t ampere_;
     GodunovFluxes_t godunov_;
     ToPrimitiveConverter_t to_primitive_;
     ToConservativeConverter_t to_conservative_;
-    FiniteVolumeEuler_t finite_volume_euler_;
-    ConstrainedTransport_t constrained_transport_;
-    Faraday_t faraday_;
+    TimeIntegrator_t time_integrator_;
 
     std::string const integrator_;
 
@@ -129,48 +125,14 @@ class SolverMHD : public ISolver<AMR_Types>
     void godunov_fluxes_(level_t& level, ModelViews_t& views, Messenger& fromCoarser,
                          double const currentTime, double const newTime);
 
-    void time_integrator_(level_t& level, ModelViews_t& views, Messenger& fromCoarser,
-                          double const currentTime, double const newTime);
+    void time_integration_(level_t& level, ModelViews_t& views, Messenger& fromCoarser,
+                           double const currentTime, double const newTime);
 
     template<typename Layout, typename View, typename VecField, typename Qs, typename... Fluxes>
     void euler_(Messenger& fromCoarser, level_t& level, double const newTime, Layout& layouts,
                 Qs quantities, View E, VecField Emodel, Qs quantities_new, double const dt,
                 Fluxes... fluxes);
 
-    template<typename Field, typename VecField>
-    struct Q
-    {
-        template<auto direction>
-        auto static rhoV(VecField const& rhoVvec)
-        {
-            Field rhoV_comp;
-            for (size_t i = 0; i < rhoVvec.size(); ++i)
-                rhoV_comp.push_back(&((*rhoVvec[i])(direction)));
-            return rhoV_comp;
-        }
-
-        Q(Field& rho_, VecField& rhoV_, VecField& B_, Field& Etot_)
-            : rho(rho_)
-            , B(B_)
-            , Etot(Etot_)
-            , rhoVx{rhoV<core::Component::X>(rhoV_)}
-            , rhoVy{rhoV<core::Component::Y>(rhoV_)}
-            , rhoVz{rhoV<core::Component::Z>(rhoV_)}
-        {
-        }
-
-        Field& rho;
-        VecField& B;
-        Field& Etot;
-        Field rhoVx, rhoVy, rhoVz;
-    };
-
-    template<typename Field, typename VecField>
-    auto static make_q(Field& rho, VecField& rhoV, VecField& B, Field& Etot)
-    {
-        return Q<Field, VecField>{rho, rhoV, B, Etot};
-    }
-
     struct TimeSetter
     {
         /*template <typename QuantityAccessor>*/
@@ -199,7 +161,7 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::advanceLevel(
 
     godunov_fluxes_(*level, modelView, fromCoarser, currentTime, newTime);
 
-    time_integrator_(*level, modelView, fromCoarser, currentTime, newTime);
+    time_integration_(*level, modelView, fromCoarser, currentTime, newTime);
 }
 
 
@@ -243,7 +205,7 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::godunov_fluxes_(
 
 
 template<typename MHDModel, typename AMR_Types, typename Messenger, typename ModelViews_t>
-void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::time_integrator_(
+void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::time_integration_(
     level_t& level, ModelViews_t& views, Messenger& fromCoarser, double const currentTime,
     double const newTime)
 {
@@ -253,118 +215,74 @@ void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::time_integrator_(
 
     to_conservative_(views.layouts, views.rho, views.V, views.B, views.P, views.rhoV, views.Etot);
 
-    Q Un = make_q(views.rho, views.rhoV, views.B, views.Etot);
-
-    Q F_x = make_q(views.rho_x, views.rhoV_x, views.B_x, views.Etot_x);
+    fromCoarser.fillMagneticFluxGhosts(views.B_x, level, newTime);
 
     if constexpr (dimension == 1)
     {
         if (integrator_ == "euler")
-            euler_(fromCoarser, level, newTime, views.layouts, Un, views.E, views.model().state.E,
-                   Un, dt, F_x);
-        /*else if (integrator_ == "TVDRK2")*/
-        /*{*/
-        /*    Q U1(views.rho1, views.rhoV1, views.B1, views.Etot1);*/
-        /*    euler_(fromCoarser, level, newTime, views.layouts, Un, views.E,
-         * views.model().state.E,*/
-        /*           U1, dt, F_x);*/
-        /*}*/
+            time_integrator_.euler(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                   views.E, dt, views.rho_x, views.rhoV_x, views.B_x, views.Etot_x);
+        else if (integrator_ == "tvdrk2")
+            time_integrator_.tvdrk2(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                    views.rho1, views.rhoV1, views.B1, views.Etot1, views.E, dt,
+                                    views.rho_x, views.rhoV_x, views.B_x, views.Etot_x);
+        else if (integrator_ == "tvdrk3")
+            time_integrator_.tvdrk3(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                    views.rho1, views.rhoV1, views.B1, views.Etot1, views.rho2,
+                                    views.rhoV2, views.B2, views.Etot2, views.E, dt, views.rho_x,
+                                    views.rhoV_x, views.B_x, views.Etot_x);
     }
     if constexpr (dimension >= 2)
     {
-        Q F_y = make_q(views.rho_y, views.rhoV_y, views.B_y, views.Etot_y);
+        fromCoarser.fillMagneticFluxGhosts(views.B_y, level, newTime);
 
         if constexpr (dimension == 2)
         {
-            euler_(fromCoarser, level, newTime, views.layouts, Un, views.E, views.model().state.E,
-                   Un, dt, F_x, F_y);
+            if (integrator_ == "euler")
+                time_integrator_.euler(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                       views.E, dt, views.rho_x, views.rhoV_x, views.B_x,
+                                       views.Etot_x, views.rho_y, views.rhoV_y, views.B_y,
+                                       views.Etot_y);
+            else if (integrator_ == "tvdrk2")
+                time_integrator_.tvdrk2(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                        views.rho1, views.rhoV1, views.B1, views.Etot1, views.E, dt,
+                                        views.rho_x, views.rhoV_x, views.B_x, views.Etot_x,
+                                        views.rho_y, views.rhoV_y, views.B_y, views.Etot_y);
+            else if (integrator_ == "tvdrk3")
+                time_integrator_.tvdrk3(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                        views.rho1, views.rhoV1, views.B1, views.Etot1, views.rho2,
+                                        views.rhoV2, views.B2, views.Etot2, views.E, dt,
+                                        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)
         {
-            Q F_z = make_q(views.rho_z, views.rhoV_z, views.B_z, views.Etot_z);
-
-            euler_(fromCoarser, level, newTime, views.layouts, Un, views.E, views.model().state.E,
-                   Un, dt, F_x, F_y, F_z);
+            fromCoarser.fillMagneticFluxGhosts(views.B_z, level, newTime);
+
+            if (integrator_ == "euler")
+                time_integrator_.euler(
+                    views.layouts, views.rho, views.rhoV, views.B, views.Etot, views.E, dt,
+                    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);
+            else if (integrator_ == "tvdrk2")
+                time_integrator_.tvdrk2(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                        views.rho1, views.rhoV1, views.B1, views.Etot1, views.E, dt,
+                                        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);
+            else if (integrator_ == "tvdrk3")
+                time_integrator_.tvdrk3(views.layouts, views.rho, views.rhoV, views.B, views.Etot,
+                                        views.rho1, views.rhoV1, views.B1, views.Etot1, views.rho2,
+                                        views.rhoV2, views.B2, views.Etot2, views.E, dt,
+                                        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);
         }
     }
 
     to_primitive_(views.layouts, views.rho, views.rhoV, views.B, views.Etot, views.V, views.P);
 }
 
-template<typename MHDModel, typename AMR_Types, typename Messenger, typename ModelViews_t>
-template<typename Layout, typename View, typename VecField, typename Qs, typename... Fluxes>
-void SolverMHD<MHDModel, AMR_Types, Messenger, ModelViews_t>::euler_(
-    Messenger& fromCoarser, level_t& level, double const newTime, Layout& layouts, Qs quantities,
-    View E, VecField Emodel, Qs quantities_new, double const dt, Fluxes... fluxes)
-{
-    auto&& flux_tuple = std::forward_as_tuple(fluxes...);
-
-    if constexpr (dimension == 1)
-    {
-        auto&& fluxes_x = std::get<0>(flux_tuple);
-
-        finite_volume_euler_(layouts, quantities.rho, quantities_new.rho, dt, fluxes_x.rho);
-        finite_volume_euler_(layouts, quantities.rhoVx, quantities_new.rhoVx, dt, fluxes_x.rhoVx);
-        finite_volume_euler_(layouts, quantities.rhoVy, quantities_new.rhoVy, dt, fluxes_x.rhoVy);
-        finite_volume_euler_(layouts, quantities.rhoVz, quantities_new.rhoVz, dt, fluxes_x.rhoVz);
-        finite_volume_euler_(layouts, quantities.Etot, quantities_new.Etot, dt, fluxes_x.Etot);
-
-        fromCoarser.fillMagneticFluxGhosts(fluxes_x.B, level, newTime);
-        constrained_transport_(layouts, E, fluxes_x.B);
-        fromCoarser.fillElectricGhosts(Emodel, level, newTime);
-
-        faraday_(layouts, quantities.B, E, quantities_new.B, dt);
-    }
-    if constexpr (dimension == 2)
-    {
-        auto&& fluxes_x = std::get<0>(flux_tuple);
-        auto&& fluxes_y = std::get<1>(flux_tuple);
-
-        finite_volume_euler_(layouts, quantities.rho, quantities_new.rho, dt, fluxes_x.rho,
-                             fluxes_y.rho);
-        finite_volume_euler_(layouts, quantities.rhoVx, quantities_new.rhoVx, dt, fluxes_x.rhoVx,
-                             fluxes_y.rhoVx);
-        finite_volume_euler_(layouts, quantities.rhoVy, quantities_new.rhoVy, dt, fluxes_x.rhoVy,
-                             fluxes_y.rhoVy);
-        finite_volume_euler_(layouts, quantities.rhoVz, quantities_new.rhoVz, dt, fluxes_x.rhoVz,
-                             fluxes_y.rhoVz);
-        finite_volume_euler_(layouts, quantities.Etot, quantities_new.Etot, dt, fluxes_x.Etot,
-                             fluxes_y.Etot);
-
-        fromCoarser.fillMagneticFluxGhosts(fluxes_x.B, level, newTime);
-        fromCoarser.fillMagneticFluxGhosts(fluxes_y.B, level, newTime);
-        constrained_transport_(layouts, E, fluxes_x.B, fluxes_y.B);
-        fromCoarser.fillElectricGhosts(Emodel, level, newTime);
-
-        faraday_(layouts, quantities.B, E, quantities_new.B, dt);
-    }
-    if constexpr (dimension == 3)
-    {
-        auto&& fluxes_x = std::get<0>(flux_tuple);
-        auto&& fluxes_y = std::get<1>(flux_tuple);
-        auto&& fluxes_z = std::get<2>(flux_tuple);
-
-        finite_volume_euler_(layouts, quantities.rho, quantities_new.rho, dt, fluxes_x.rho,
-                             fluxes_y.rho, fluxes_z.rho);
-        finite_volume_euler_(layouts, quantities.rhoVx, quantities_new.rhoVx, dt, fluxes_x.rhoVx,
-                             fluxes_y.rhoVx, fluxes_z.rhoVx);
-        finite_volume_euler_(layouts, quantities.rhoVy, quantities_new.rhoVy, dt, fluxes_x.rhoVy,
-                             fluxes_y.rhoVy, fluxes_z.rhoVy);
-        finite_volume_euler_(layouts, quantities.rhoVz, quantities_new.rhoVz, dt, fluxes_x.rhoVz,
-                             fluxes_y.rhoVz, fluxes_z.rhoVz);
-        finite_volume_euler_(layouts, quantities.Etot, quantities_new.Etot, dt, fluxes_x.Etot,
-                             fluxes_y.Etot, fluxes_z.Etot);
-
-        fromCoarser.fillMagneticFluxGhosts(fluxes_x.B, level, newTime);
-        fromCoarser.fillMagneticFluxGhosts(fluxes_y.B, level, newTime);
-        fromCoarser.fillMagneticFluxGhosts(fluxes_z.B, level, newTime);
-        constrained_transport_(layouts, E, fluxes_x.B, fluxes_y.B, fluxes_z.B);
-        fromCoarser.fillElectricGhosts(Emodel, level, newTime);
-
-        faraday_(layouts, quantities.B, E, quantities_new.B, dt);
-    }
-}
-
 } // namespace PHARE::solver
 
 #endif

src/amr/solvers/solver_mhd_model_view.hpp

@@ -9,6 +9,7 @@
 #include "core/numerics/finite_volume_euler/finite_volume_euler.hpp"
 #include "core/numerics/primite_conservative_converter/to_conservative_converter.hpp"
 #include "core/numerics/primite_conservative_converter/to_primitive_converter.hpp"
+#include "core/numerics/time_integrator/time_integrator.hpp"
 
 namespace PHARE::solver
 {
@@ -118,6 +119,57 @@ class ToPrimitiveTransformer
     core_type to_primtitve_;
 };
 
+template<typename GridLayout>
+class TimeIntegratorTransformer
+{
+    using core_type = PHARE::core::TimeIntegrator<GridLayout>;
+
+public:
+    template<typename Layout, typename Field, typename VecField, typename... Fluxes>
+    void euler(Layout const& layouts, Field& rho, VecField& rhoV, VecField& B, Field& Etot,
+               VecField& E, double const dt, Fluxes&... fluxes)
+    {
+        assert_equal_sizes(rho, rhoV, B, Etot, E, fluxes...);
+        for (std::size_t i = 0; i < layouts.size(); ++i)
+        {
+            auto _ = core::SetLayout(layouts[i], time_integrator_);
+            time_integrator_.euler(*rho[i], *rhoV[i], *B[i], *Etot[i], *E[i], dt, *fluxes[i]...);
+        }
+    }
+
+    template<typename Layout, typename Field, typename VecField, typename... Fluxes>
+    void tvdrk2(Layout const& layouts, Field& rho, VecField& rhoV, VecField& B, Field& Etot,
+                Field& rho1, VecField& rhoV1, VecField& B1, Field& Etot1, VecField& E,
+                double const dt, Fluxes&... fluxes)
+    {
+        assert_equal_sizes(rho, rhoV, B, Etot, E, fluxes...);
+        for (std::size_t i = 0; i < layouts.size(); ++i)
+        {
+            auto _ = core::SetLayout(layouts[i], time_integrator_);
+            time_integrator_.tvdrk2(*rho[i], *rhoV[i], *B[i], *Etot[i], *rho1[i], *rhoV1[i], *B1[i],
+                                    *Etot1[i], *E[i], dt, *fluxes[i]...);
+        }
+    }
+
+    template<typename Layout, typename Field, typename VecField, typename... Fluxes>
+    void tvdrk3(Layout const& layouts, Field& rho, VecField& rhoV, VecField& B, Field& Etot,
+                Field& rho1, VecField& rhoV1, VecField& B1, Field& Etot1, Field& rho2,
+                VecField& rhoV2, VecField& B2, Field& Etot2, VecField& E, double const dt,
+                Fluxes&... fluxes)
+    {
+        assert_equal_sizes(rho, rhoV, B, Etot, E, fluxes...);
+        for (std::size_t i = 0; i < layouts.size(); ++i)
+        {
+            auto _ = core::SetLayout(layouts[i], time_integrator_);
+            time_integrator_.tvdrk3(*rho[i], *rhoV[i], *B[i], *Etot[i], *rho1[i], *rhoV1[i], *B1[i],
+                                    *Etot1[i], *rho2[i], *rhoV2[i], *B2[i], *Etot2[i], *E[i], dt,
+                                    *fluxes[i]...);
+        }
+    }
+
+    core_type time_integrator_;
+};
+
 
 template<typename MHDModel_>
 class MHDModelView : public ISolverModelView
@@ -127,11 +179,13 @@ class MHDModelView : public ISolverModelView
     using GridLayout                = typename MHDModel_t::gridlayout_type;
     using GodunovFluxes_t           = GodunovFluxesTransformer<GridLayout>;
     using Ampere_t                  = AmpereTransformer<GridLayout>;
-    using FiniteVolumeEuler_t       = FiniteVolumeEulerTransformer<GridLayout>;
-    using ConstrainedTransport_t    = ConstrainedTransportTransformer<GridLayout>;
-    using Faraday_t                 = FaradayTransformer<GridLayout>;
     using ToPrimitiveConverter_t    = ToPrimitiveTransformer<GridLayout>;
     using ToConservativeConverter_t = ToConservativeTransformer<GridLayout>;
+    using TimeIntegrator_t          = TimeIntegratorTransformer<GridLayout>;
+
+    using FiniteVolumeEuler_t    = FiniteVolumeEulerTransformer<GridLayout>;
+    using ConstrainedTransport_t = ConstrainedTransportTransformer<GridLayout>;
+    using Faraday_t              = FaradayTransformer<GridLayout>;
 };
 
 }; // namespace PHARE::solver

src/core/data/grid/gridlayout.hpp

@@ -566,7 +566,7 @@ namespace core
          * on the dimensionality of the GridLayout.
          */
         template<auto direction, typename Field>
-        NO_DISCARD auto deriv(Field const& operand, MeshIndex<Field::dimension> index)
+        NO_DISCARD auto deriv(Field const& operand, MeshIndex<Field::dimension> index) const
         {
             auto fieldCentering = centering(operand.physicalQuantity());
             using PHARE::core::dirX;