Skip to content

Commit 4239b27

Browse files
UCaromelUCaromel
committed
implemented SSPRK4_5
1 parent 813b83e commit 4239b27

File tree

8 files changed

+341
-97
lines changed

8 files changed

+341
-97
lines changed

‎src/amr/solvers/solver_mhd.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -420,8 +420,8 @@ void SolverMHD<MHDModel, AMR_Types, TimeIntegratorStrategy, Messenger, ModelView
420420
auto& mhdModel = dynamic_cast<MHDModel&>(model);
421421
auto&& [timeFluxes, timeElectric] = evolve_.exposeFluxes();
422422

423-
reflux_euler_(mhdModel, stateOld_, mhdModel.state, timeElectric, timeFluxes, bc, level,
424-
oldTime_[level.getLevelNumber()], time);
423+
reflux_euler_(mhdModel, stateOld_, mhdModel.state, timeElectric, timeFluxes, bc, level, time,
424+
time - oldTime_[level.getLevelNumber()]);
425425
}
426426

427427
template<typename MHDModel, typename AMR_Types, typename TimeIntegratorStrategy, typename Messenger,

‎src/amr/solvers/time_integrator/compute_fluxes.hpp‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
#ifndef PHARE_CORE_NUMERICS_TIME_INTEGRATOR_COMPUTE_FLUXES_HPP
2+
#define PHARE_CORE_NUMERICS_TIME_INTEGRATOR_COMPUTE_FLUXES_HPP
3+
4+
#include "initializer/data_provider.hpp"
5+
#include "amr/solvers/solver_mhd_model_view.hpp"
6+
7+
namespace PHARE::solver
8+
{
9+
template<template<typename> typename FVMethodStrategy, typename MHDModel>
10+
class ComputeFluxes
11+
{
12+
using level_t = typename MHDModel::level_t;
13+
using Layout = typename MHDModel::gridlayout_type;
14+
using Dispatchers_t = Dispatchers<Layout>;
15+
16+
using Ampere_t = Dispatchers_t::Ampere_t;
17+
using FVMethod_t = Dispatchers_t::template FVMethod_t<FVMethodStrategy>;
18+
using ConstrainedTransport_t = Dispatchers_t::ConstrainedTransport_t;
19+
20+
using ToPrimitiveConverter_t = Dispatchers_t::ToPrimitiveConverter_t;
21+
using ToConservativeConverter_t = Dispatchers_t::ToConservativeConverter_t;
22+
23+
constexpr static auto Hall = FVMethod_t::Hall;
24+
constexpr static auto Resistivity = FVMethod_t::Resistivity;
25+
constexpr static auto HyperResistivity = FVMethod_t::HyperResistivity;
26+
27+
public:
28+
ComputeFluxes(PHARE::initializer::PHAREDict const& dict)
29+
: fvm_{dict["fv_method"]}
30+
, to_primitive_{dict["to_primitive"]}
31+
, to_conservative_{dict["to_conservative"]}
32+
{
33+
}
34+
35+
void operator()(MHDModel& model, auto& state, auto& fluxes, auto& bc, level_t& level,
36+
double const newTime)
37+
{
38+
to_primitive_(level, model, newTime, state);
39+
40+
if constexpr (Hall || Resistivity || HyperResistivity)
41+
{
42+
ampere_(level, model, newTime, state);
43+
44+
bc.fillCurrentGhosts(state.J, level, newTime);
45+
}
46+
47+
fvm_(level, model, newTime, state, fluxes);
48+
49+
// unecessary if we decide to store both primitive and conservative variables
50+
to_conservative_(level, model, newTime, state);
51+
52+
bc.fillMagneticFluxesXGhosts(fluxes.B_fx, level, newTime);
53+
54+
if constexpr (MHDModel::dimension >= 2)
55+
{
56+
bc.fillMagneticFluxesYGhosts(fluxes.B_fy, level, newTime);
57+
58+
if constexpr (MHDModel::dimension == 3)
59+
{
60+
bc.fillMagneticFluxesZGhosts(fluxes.B_fz, level, newTime);
61+
}
62+
}
63+
64+
ct_(level, model, state, fluxes);
65+
66+
bc.fillElectricGhosts(state.E, level, newTime);
67+
}
68+
69+
private:
70+
Ampere_t ampere_;
71+
FVMethod_t fvm_;
72+
ConstrainedTransport_t ct_;
73+
ToPrimitiveConverter_t to_primitive_;
74+
ToConservativeConverter_t to_conservative_;
75+
};
76+
} // namespace PHARE::solver
77+
78+
#endif

‎src/amr/solvers/time_integrator/euler.hpp‎

Lines changed: 14 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -2,88 +2,40 @@
22
#define PHARE_CORE_NUMERICS_TIME_INTEGRATOR_EULER_HPP
33

44
#include "initializer/data_provider.hpp"
5-
#include "amr/solvers/solver_mhd_model_view.hpp"
5+
#include "amr/solvers/time_integrator/compute_fluxes.hpp"
6+
#include "amr/solvers/time_integrator/euler_using_computed_flux.hpp"
67

78
namespace PHARE::solver
89
{
910
template<template<typename> typename FVMethodStrategy, typename MHDModel>
1011
class Euler
1112
{
12-
using level_t = typename MHDModel::level_t;
13-
using Layout = typename MHDModel::gridlayout_type;
14-
using Dispatchers_t = Dispatchers<Layout>;
13+
using level_t = typename MHDModel::level_t;
1514

16-
using Ampere_t = Dispatchers_t::Ampere_t;
17-
using FVMethod_t = Dispatchers_t::template FVMethod_t<FVMethodStrategy>;
18-
using FiniteVolumeEuler_t = Dispatchers_t::FiniteVolumeEuler_t;
19-
using ConstrainedTransport_t = Dispatchers_t::ConstrainedTransport_t;
20-
using Faraday_t = Dispatchers_t::Faraday_t;
21-
22-
using ToPrimitiveConverter_t = Dispatchers_t::ToPrimitiveConverter_t;
23-
using ToConservativeConverter_t = Dispatchers_t::ToConservativeConverter_t;
24-
25-
constexpr static auto Hall = FVMethod_t::Hall;
26-
constexpr static auto Resistivity = FVMethod_t::Resistivity;
27-
constexpr static auto HyperResistivity = FVMethod_t::HyperResistivity;
15+
using ComputeFluxes_t = ComputeFluxes<FVMethodStrategy, MHDModel>;
16+
using EulerUsingComputedFlux_t = EulerUsingComputedFlux<MHDModel>;
2817

2918
public:
3019
Euler(PHARE::initializer::PHAREDict const& dict)
31-
: fvm_{dict["fv_method"]}
32-
, to_primitive_{dict["to_primitive"]}
33-
, to_conservative_{dict["to_conservative"]}
20+
: compute_fluxes_{dict}
3421
{
3522
}
3623

3724
void operator()(MHDModel& model, auto& state, auto& statenew, auto& fluxes, auto& bc,
38-
level_t& level, double const currentTime, double const newTime)
25+
level_t& level, double const currentTime, double const newTime,
26+
double dt = std::nan(""))
3927
{
40-
double const dt = newTime - currentTime;
41-
42-
to_primitive_(level, model, newTime, state);
43-
44-
if constexpr (Hall || Resistivity || HyperResistivity)
45-
{
46-
ampere_(level, model, newTime, state);
47-
48-
bc.fillCurrentGhosts(state.J, level, newTime);
49-
}
50-
51-
fvm_(level, model, newTime, state, fluxes);
52-
53-
// unecessary if we decide to store both primitive and conservative variables
54-
to_conservative_(level, model, newTime, state);
55-
56-
fv_euler_(level, model, newTime, state, statenew, fluxes, dt);
57-
58-
bc.fillMagneticFluxesXGhosts(fluxes.B_fx, level, newTime);
59-
60-
if constexpr (MHDModel::dimension >= 2)
61-
{
62-
bc.fillMagneticFluxesYGhosts(fluxes.B_fy, level, newTime);
63-
64-
if constexpr (MHDModel::dimension == 3)
65-
{
66-
bc.fillMagneticFluxesZGhosts(fluxes.B_fz, level, newTime);
67-
}
68-
}
69-
70-
ct_(level, model, state, fluxes);
71-
72-
bc.fillElectricGhosts(state.E, level, newTime);
28+
if (std::isnan(dt))
29+
dt = newTime - currentTime;
7330

74-
faraday_(level, model, state, state.E, statenew, dt);
31+
compute_fluxes_(model, state, fluxes, bc, level, newTime);
7532

76-
bc.fillMomentsGhosts(statenew, level, newTime);
33+
euler_using_computed_flux_(model, state, statenew, state.E, fluxes, bc, level, newTime, dt);
7734
}
7835

7936
private:
80-
Ampere_t ampere_;
81-
FVMethod_t fvm_;
82-
FiniteVolumeEuler_t fv_euler_;
83-
ConstrainedTransport_t ct_;
84-
Faraday_t faraday_;
85-
ToPrimitiveConverter_t to_primitive_;
86-
ToConservativeConverter_t to_conservative_;
37+
ComputeFluxes_t compute_fluxes_;
38+
EulerUsingComputedFlux_t euler_using_computed_flux_;
8739
};
8840
} // namespace PHARE::solver
8941

‎src/amr/solvers/time_integrator/euler_using_computed_flux.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,11 +19,11 @@ class EulerUsingComputedFlux
1919
public:
2020
EulerUsingComputedFlux() {}
2121

22+
// we provide dt here because we sometimes need it to be different from newTime-currentTime, for
23+
// example in the case of some rk integration methods
2224
void operator()(MHDModel& model, auto& state, auto& statenew, auto& E, auto& fluxes, auto& bc,
23-
level_t& level, double const currentTime, double const newTime)
25+
level_t& level, double const newTime, double const dt)
2426
{
25-
double const dt = newTime - currentTime;
26-
2727
fv_euler_(level, model, newTime, state, statenew, fluxes, dt);
2828

2929
faraday_(level, model, state, E, statenew, dt);

‎src/amr/solvers/time_integrator/ssprk4_5_integrator.hpp‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
#ifndef PHARE_CORE_NUMERICS_SSPRK4_5_INTEGRATOR_HPP
2+
#define PHARE_CORE_NUMERICS_SSPRK4_5_INTEGRATOR_HPP
3+
4+
#include "initializer/data_provider.hpp"
5+
#include "amr/solvers/time_integrator/base_mhd_timestepper.hpp"
6+
#include "amr/solvers/time_integrator/compute_fluxes.hpp"
7+
#include "amr/solvers/time_integrator/euler_using_computed_flux.hpp"
8+
#include "amr/solvers/solver_mhd_model_view.hpp"
9+
#include "amr/solvers/time_integrator/euler.hpp"
10+
#include "core/numerics/time_integrator_utils.hpp"
11+
12+
namespace PHARE::solver
13+
{
14+
template<template<typename> typename FVMethodStrategy, typename MHDModel>
15+
class SSPRK4_5Integrator : public BaseMHDTimestepper<MHDModel>
16+
{
17+
using Super = BaseMHDTimestepper<MHDModel>;
18+
19+
using level_t = typename MHDModel::level_t;
20+
using FieldT = typename MHDModel::field_type;
21+
using VecFieldT = typename MHDModel::vecfield_type;
22+
using GridLayoutT = typename MHDModel::gridlayout_type;
23+
using MHDStateT = typename MHDModel::state_type;
24+
25+
using Dispatchers_t = Dispatchers<GridLayoutT>;
26+
using RKUtils_t = Dispatchers_t::RKUtils_t;
27+
28+
using RKPair_t = core::RKPair<typename VecFieldT::value_type, MHDStateT>;
29+
30+
public:
31+
SSPRK4_5Integrator(PHARE::initializer::PHAREDict const& dict)
32+
: Super{dict}
33+
, euler_{dict}
34+
, compute_fluxes_{dict}
35+
{
36+
}
37+
38+
// Butcher fluxes are used to accumulate fluxes over multiple stages, the corresponding buffer
39+
// should only contain the fluxes over one time step. The accumulation over all substeps is
40+
// delegated to the solver.
41+
void operator()(MHDModel& model, auto& state, auto& fluxes, auto& bc, level_t& level,
42+
double const currentTime, double const newTime)
43+
{
44+
this->resetButcherFluxes_(model, level);
45+
46+
auto const dt = newTime - currentTime;
47+
48+
// U1 = Un + w0_*dt*F(Un)
49+
euler_(model, state, state1_, fluxes, bc, level, currentTime, newTime, w0_ * dt);
50+
51+
this->accumulateButcherFluxes_(
52+
model, state.E, fluxes, level,
53+
(w0_ * w11_ * w21_ * w31_ * w43_ + w0_ * w11_ * w21_ * w41_ + w0_ * w11_));
54+
55+
// U2 = w10_*Un + w11_*U1 + w12_*dt*F(U1)
56+
//
57+
// U2 = w10_Un + w11_*U1
58+
rk_step_(level, model, newTime, state2_, RKPair_t{w10_, state}, RKPair_t{w11_, state1_});
59+
60+
// U2 = U2 + w12_*dt*F(U1)
61+
compute_fluxes_(model, state1_, fluxes, bc, level, newTime);
62+
63+
euler_using_butcher_fluxes_(model, state2_, state2_, state1_.E, fluxes, bc, level, newTime,
64+
w12_ * dt);
65+
66+
this->accumulateButcherFluxes_(
67+
model, state1_.E, fluxes, level,
68+
(w12_ * w21_ * w31_ * w43_ + w12_ * w21_ * w41_ + w12_ * w40_));
69+
70+
// U3 = w20_*Un + w21_*U2 + w22_*dt*F(U2)
71+
//
72+
// U3 = w20_*Un + w21_*U2
73+
rk_step_(level, model, newTime, state3_, RKPair_t{w20_, state}, RKPair_t{w21_, state2_});
74+
75+
// U3 = U3 + w22_*dt*F(U2)
76+
compute_fluxes_(model, state2_, fluxes, bc, level, newTime);
77+
78+
euler_using_butcher_fluxes_(model, state3_, state3_, state2_.E, fluxes, bc, level, newTime,
79+
w22_ * dt);
80+
81+
this->accumulateButcherFluxes_(model, state2_.E, fluxes, level,
82+
(w22_ * w31_ * w43_ + w22_ * w41_));
83+
84+
// U4 = w30_*Un + w31_*U3 + w32_*dt*F(U3)
85+
//
86+
// U4 = w30_*Un + w31_*U3
87+
rk_step_(level, model, newTime, state4_, RKPair_t{w30_, state}, RKPair_t{w31_, state3_});
88+
89+
// U4 = U4 + w32_*dt*F(U3)
90+
// if we were not using butcher formulation, we would need a separate flux buffer for F(U3)
91+
// for the final step
92+
compute_fluxes_(model, state3_, fluxes, bc, level, newTime);
93+
94+
euler_using_butcher_fluxes_(model, state4_, state4_, state3_.E, fluxes, bc, level, newTime,
95+
w32_ * dt);
96+
97+
this->accumulateButcherFluxes_(model, state3_.E, fluxes, level, (w32_ * w43_ + w42_));
98+
99+
compute_fluxes_(model, state4_, fluxes, bc, level, newTime);
100+
101+
this->accumulateButcherFluxes_(model, state4_.E, fluxes, level, w44_);
102+
103+
euler_using_butcher_fluxes_(model, state, state, this->butcherE_, this->butcherFluxes_, bc,
104+
level, newTime, dt);
105+
106+
// Un+1 = w40_*U2 + w41_*U3 + w42_*F(U3) + w43_*U4 + w44_*dt*F(U4)
107+
}
108+
109+
void registerResources(MHDModel& model)
110+
{
111+
Super::registerResources(model);
112+
model.resourcesManager->registerResources(state1_);
113+
model.resourcesManager->registerResources(state2_);
114+
}
115+
116+
void allocate(MHDModel& model, auto& patch, double const allocateTime) const
117+
{
118+
Super::allocate(model, patch, allocateTime);
119+
model.resourcesManager->allocate(state1_, patch, allocateTime);
120+
model.resourcesManager->allocate(state2_, patch, allocateTime);
121+
}
122+
123+
void fillMessengerInfo(auto& info) const
124+
{
125+
info.ghostDensity.push_back(state1_.rho.name());
126+
info.ghostMomentum.push_back(state1_.rhoV.name());
127+
info.ghostTotalEnergy.push_back(state1_.Etot.name());
128+
info.ghostElectric.push_back(state1_.E.name());
129+
info.ghostCurrent.push_back(state1_.J.name());
130+
131+
info.ghostDensity.push_back(state2_.rho.name());
132+
info.ghostMomentum.push_back(state2_.rhoV.name());
133+
info.ghostTotalEnergy.push_back(state2_.Etot.name());
134+
info.ghostElectric.push_back(state2_.E.name());
135+
info.ghostCurrent.push_back(state2_.J.name());
136+
}
137+
138+
NO_DISCARD auto getCompileTimeResourcesViewList()
139+
{
140+
return std::tuple_cat(Super::getCompileTimeResourcesViewList(),
141+
std::forward_as_tuple(state1_, state2_));
142+
}
143+
144+
NO_DISCARD auto getCompileTimeResourcesViewList() const
145+
{
146+
return std::tuple_cat(Super::getCompileTimeResourcesViewList(),
147+
std::forward_as_tuple(state1_, state2_));
148+
}
149+
150+
using Super::exposeFluxes;
151+
152+
private:
153+
static constexpr auto w0_{0.391752226571890};
154+
static constexpr auto w10_{0.444370493651235};
155+
static constexpr auto w11_{0.555629506348765};
156+
static constexpr auto w12_{0.368410593050371};
157+
static constexpr auto w20_{0.620101851488403};
158+
static constexpr auto w21_{0.379898148511597};
159+
static constexpr auto w22_{0.251891774271694};
160+
static constexpr auto w30_{0.178079954393132};
161+
static constexpr auto w31_{0.821920045606868};
162+
static constexpr auto w32_{0.544974750228521};
163+
static constexpr auto w40_{0.517231671970585};
164+
static constexpr auto w41_{0.096059710526147};
165+
static constexpr auto w42_{0.063692468666290};
166+
static constexpr auto w43_{0.386708617503268};
167+
static constexpr auto w44_{0.226007483236906};
168+
169+
Euler<FVMethodStrategy, MHDModel> euler_;
170+
ComputeFluxes<FVMethodStrategy, MHDModel> compute_fluxes_;
171+
EulerUsingComputedFlux<MHDModel> euler_using_butcher_fluxes_;
172+
RKUtils_t rk_step_;
173+
174+
MHDStateT state1_{"state1"};
175+
MHDStateT state2_{"state2"};
176+
MHDStateT state3_{"state2"};
177+
MHDStateT state4_{"state2"};
178+
};
179+
180+
} // namespace PHARE::solver
181+
182+
#endif

‎src/amr/solvers/time_integrator/tvdrk2_integrator.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ class TVDRK2Integrator : public BaseMHDTimestepper<MHDModel>
5151
this->accumulateButcherFluxes_(model, state1_.E, fluxes, level, w1_);
5252

5353
euler_using_butcher_fluxes_(model, state, state, this->butcherE_, this->butcherFluxes_, bc,
54-
level, currentTime, newTime);
54+
level, newTime, newTime - currentTime);
5555

5656
// Un+1 = 0.5*Un + 0.5*Euler(U1)
5757
// tvdrk2_step_(level, model, newTime, state, RKPair_t{w0_, state}, RKPair_t{w1_, state1_});

0 commit comments

Comments
 (0)