Skip to content

Commit 05eac19

Browse files
committed
fixed added setTime for time refiners
1 parent ba84522 commit 05eac19

File tree

14 files changed

+499
-156
lines changed

14 files changed

+499
-156
lines changed

pyphare/pyphare/pharein/simulation.py

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -652,6 +652,13 @@ def check_mhd_constants(**kwargs):
652652

653653
return gamma, eta, nu
654654

655+
def check_mhd_terms(**kwargs):
656+
hall = kwargs.get("hall", False)
657+
res = kwargs.get("res", False)
658+
hyper_res = kwargs.get("hyper_res", False)
659+
660+
return hall, res, hyper_res
661+
655662
def check_mhd_parameters(**kwargs):
656663
reconstruction = kwargs.get("reconstruction", "")
657664
limiter = kwargs.get("limiter", "")
@@ -704,6 +711,9 @@ def wrapper(simulation_object, **kwargs):
704711
"gamma",
705712
"eta",
706713
"nu",
714+
"hall",
715+
"res",
716+
"hyper_res",
707717
"reconstruction",
708718
"limiter",
709719
"riemann",
@@ -793,6 +803,11 @@ def wrapper(simulation_object, **kwargs):
793803
kwargs["eta"] = eta
794804
kwargs["nu"] = nu
795805

806+
hall, res, hyper_res = check_mhd_terms(**kwargs)
807+
kwargs["hall"] = hall
808+
kwargs["res"] = res
809+
kwargs["hyper_res"] = hyper_res
810+
796811
reconstruction, limiter, riemann, mhd_timestepper = check_mhd_parameters(**kwargs)
797812
kwargs["reconstruction"] = reconstruction
798813
kwargs["limiter"] = limiter

src/amr/messengers/mhd_messenger.hpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -139,8 +139,8 @@ namespace amr
139139
std::shared_ptr<SAMRAI::hier::PatchLevel> const& oldLevel,
140140
IPhysicalModel& model, double const initDataTime) override
141141
{
142-
auto& hybridModel = static_cast<MHDModel&>(model);
143-
auto level = hierarchy->getPatchLevel(levelNumber);
142+
auto& mhdModel = static_cast<MHDModel&>(model);
143+
auto level = hierarchy->getPatchLevel(levelNumber);
144144

145145
bool isRegriddingL0 = levelNumber == 0 and oldLevel;
146146

@@ -151,7 +151,7 @@ namespace amr
151151

152152
if (!isRegriddingL0)
153153
{
154-
auto& B = hybridModel.state.B;
154+
auto& B = mhdModel.state.B;
155155
magGhostsRefiners_.fill(B, levelNumber, initDataTime);
156156

157157
fix_magnetic_divergence_(*hierarchy, levelNumber, B);

src/amr/solvers/solver_mhd_model_view.hpp

Lines changed: 61 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -14,20 +14,42 @@
1414

1515
namespace PHARE::solver
1616
{
17+
template<typename MHDModel>
18+
struct TimeSetter
19+
{
20+
template<typename... QuantityAccessors>
21+
void operator()(auto& patch, QuantityAccessors... accessors)
22+
{
23+
(model.resourcesManager->setTime(accessors(), patch, newTime), ...);
24+
}
25+
26+
MHDModel& model;
27+
double newTime;
28+
};
29+
1730
template<typename GridLayout>
1831
class ToConservativeTransformer
1932
{
2033
using core_type = PHARE::core::ToConservativeConverter<GridLayout>;
2134

2235
public:
2336
template<typename MHDModel>
24-
void operator()(MHDModel::level_t const& level, MHDModel& model, MHDModel::state_type& state)
37+
void operator()(MHDModel::level_t const& level, MHDModel& model, double const newTime,
38+
MHDModel::state_type& state)
2539
{
40+
TimeSetter setTime{model, newTime};
41+
2642
for (auto const& patch : level)
2743
{
2844
auto layout = PHARE::amr::layoutFromPatch<GridLayout>(*patch);
2945
auto _sp = model.resourcesManager->setOnPatch(*patch, state);
3046
auto _sl = core::SetLayout(&layout, to_conservative_);
47+
48+
setTime(
49+
*patch, [&]() -> auto&& { return state.rho; }, [&]() -> auto&& { return state.V; },
50+
[&]() -> auto&& { return state.P; }, [&]() -> auto&& { return state.rhoV; },
51+
[&]() -> auto&& { return state.Etot; });
52+
3153
to_conservative_(state.rho, state.V, state.B, state.P, state.rhoV, state.Etot);
3254
}
3355
}
@@ -42,13 +64,22 @@ class ToPrimitiveTransformer
4264

4365
public:
4466
template<typename MHDModel>
45-
void operator()(MHDModel::level_t const& level, MHDModel& model, MHDModel::state_type& state)
67+
void operator()(MHDModel::level_t const& level, MHDModel& model, double const newTime,
68+
MHDModel::state_type& state)
4669
{
70+
TimeSetter setTime{model, newTime};
71+
4772
for (auto const& patch : level)
4873
{
4974
auto layout = PHARE::amr::layoutFromPatch<GridLayout>(*patch);
5075
auto _sp = model.resourcesManager->setOnPatch(*patch, state);
5176
auto _sl = core::SetLayout(&layout, to_primitive_);
77+
78+
setTime(
79+
*patch, [&]() -> auto&& { return state.rho; },
80+
[&]() -> auto&& { return state.rhoV; }, [&]() -> auto&& { return state.Etot; },
81+
[&]() -> auto&& { return state.V; }, [&]() -> auto&& { return state.P; });
82+
5283
to_primitive_(state.rho, state.rhoV, state.B, state.Etot, state.V, state.P);
5384
}
5485
}
@@ -63,14 +94,21 @@ class FVMethodTransformer
6394

6495
public:
6596
template<typename MHDModel>
66-
void operator()(MHDModel::level_t const& level, MHDModel& model, MHDModel::state_type& state,
67-
auto& fluxes)
97+
void operator()(MHDModel::level_t const& level, MHDModel& model, double const newTime,
98+
MHDModel::state_type& state, auto& fluxes)
6899
{
100+
TimeSetter setTime{model, newTime};
101+
69102
for (auto const& patch : level)
70103
{
71104
auto layout = PHARE::amr::layoutFromPatch<GridLayout>(*patch);
72105
auto _sp = model.resourcesManager->setOnPatch(*patch, state, fluxes);
73106
auto _sl = core::SetLayout(&layout, fvm_);
107+
108+
setTime(
109+
*patch, [&]() -> auto&& { return state.rho; }, [&]() -> auto&& { return state.V; },
110+
[&]() -> auto&& { return state.P; }, [&]() -> auto&& { return state.J; });
111+
74112
fvm_(state, fluxes);
75113
}
76114
}
@@ -86,14 +124,22 @@ class FiniteVolumeEulerTransformer
86124

87125
public:
88126
template<typename MHDModel>
89-
void operator()(MHDModel::level_t const& level, MHDModel& model, MHDModel::state_type& state,
90-
MHDModel::state_type& statenew, auto& fluxes, double const dt)
127+
void operator()(MHDModel::level_t const& level, MHDModel& model, double const newTime,
128+
MHDModel::state_type& state, MHDModel::state_type& statenew, auto& fluxes,
129+
double const dt)
91130
{
131+
TimeSetter setTime{model, newTime};
132+
92133
for (auto const& patch : level)
93134
{
94135
auto layout = PHARE::amr::layoutFromPatch<GridLayout>(*patch);
95136
auto _sp = model.resourcesManager->setOnPatch(*patch, state, statenew, fluxes);
96137
auto _sl = core::SetLayout(&layout, euler_);
138+
139+
setTime(
140+
*patch, [&]() -> auto&& { return state.rho; },
141+
[&]() -> auto&& { return state.rhoV; }, [&]() -> auto&& { return state.Etot; });
142+
97143
euler_(state, statenew, fluxes, dt);
98144
}
99145
}
@@ -152,14 +198,21 @@ class RKUtilsTransformer
152198

153199
public:
154200
template<typename MHDModel, typename... Pairs>
155-
void operator()(MHDModel::level_t const& level, MHDModel& model, MHDModel::state_type& res,
156-
Pairs... pairs)
201+
void operator()(MHDModel::level_t const& level, MHDModel& model, double const newTime,
202+
MHDModel::state_type& res, Pairs... pairs)
157203
{
204+
TimeSetter setTime{model, newTime};
205+
158206
for (auto const& patch : level)
159207
{
160208
auto layout = PHARE::amr::layoutFromPatch<GridLayout>(*patch);
161209
auto _sp = model.resourcesManager->setOnPatch(*patch, res, pairs.state...);
162210
auto _sl = core::SetLayout(&layout, rkutils_);
211+
212+
setTime(
213+
*patch, [&]() -> auto&& { return res.rho; }, [&]() -> auto&& { return res.rhoV; },
214+
[&]() -> auto&& { return res.Etot; });
215+
163216
rkutils_(res, pairs...);
164217
}
165218
}

src/amr/solvers/time_integrator/euler.hpp

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -35,16 +35,16 @@ class Euler
3535
{
3636
double const dt = newTime - currentTime;
3737

38-
to_primitive_(level, model, state);
38+
to_primitive_(level, model, newTime, state);
3939

4040
bc.fillMomentsGhosts(state, level.getLevelNumber(), newTime);
4141

42-
fvm_(level, model, state, fluxes);
42+
fvm_(level, model, newTime, state, fluxes);
4343

4444
// unecessary if we decide to store both primitive and conservative variables
45-
to_conservative_(level, model, state);
45+
to_conservative_(level, model, newTime, state);
4646

47-
fv_euler_(level, model, state, statenew, fluxes, dt);
47+
fv_euler_(level, model, newTime, state, statenew, fluxes, dt);
4848

4949
bc.fillMagneticFluxesXGhosts(fluxes.B_fx, level.getLevelNumber(), newTime);
5050

src/amr/solvers/time_integrator/tvdrk2_integrator.hpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ class TVDRK2Integrator
4040
euler_(model, state1_, state1_, fluxes, bc, level, currentTime, newTime);
4141

4242
// Un+1 = 0.5*Un + 0.5*Euler(U1)
43-
tvdrk2_step_(level, model, state, RKPair_t{w0_, state}, RKPair_t{w1_, state1_});
43+
tvdrk2_step_(level, model, newTime, state, RKPair_t{w0_, state}, RKPair_t{w1_, state1_});
4444
}
4545

4646
void registerResources(MHDModel& model) { model.resourcesManager->registerResources(state1_); }

src/amr/solvers/time_integrator/tvdrk3_integrator.hpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -40,13 +40,14 @@ class TVDRK3Integrator
4040
euler_(model, state1_, state1_, fluxes, bc, level, currentTime, newTime);
4141

4242
// U2 = 0.75*Un + 0.25*U1
43-
tvdrk3_step_(level, model, state2_, RKPair_t{w00_, state}, RKPair_t{w01_, state1_});
43+
tvdrk3_step_(level, model, newTime, state2_, RKPair_t{w00_, state},
44+
RKPair_t{w01_, state1_});
4445

4546
// U2 = Euler(U2)
4647
euler_(model, state2_, state2_, fluxes, bc, level, currentTime, newTime);
4748

4849
// Un+1 = 1/3*Un + 2/3*Euler(U2)
49-
tvdrk3_step_(level, model, state, RKPair_t{w10_, state}, RKPair_t{w11_, state2_});
50+
tvdrk3_step_(level, model, newTime, state, RKPair_t{w10_, state}, RKPair_t{w11_, state2_});
5051
}
5152

5253
void registerResources(MHDModel& model)

src/python3/cpp_mhd_parameters.hpp

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -217,6 +217,17 @@ constexpr void declare_all_mhd_params(py::module& m)
217217
ReconstructionType::WENO3, SlopeLimiterType::count,
218218
RiemannSolverType::Rusanov, false, false, false>::declare_etc(m, full_type);
219219

220+
variant_name = "tvdrk3_weno3_rusanov_hall";
221+
full_type = type_name + "_" + variant_name;
222+
223+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
224+
ReconstructionType::WENO3, SlopeLimiterType::count,
225+
RiemannSolverType::Rusanov, true, false, false>::declare_sim(m, full_type);
226+
227+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
228+
ReconstructionType::WENO3, SlopeLimiterType::count,
229+
RiemannSolverType::Rusanov, true, false, false>::declare_etc(m, full_type);
230+
220231
/*auto constexpr ti_tuple = make_enum_tuple<TimeIntegratorType>();*/
221232
/*auto constexpr rc_tuple = make_enum_tuple<ReconstructionType>();*/
222233
/*auto constexpr sl_tuple = make_enum_tuple<SlopeLimiterType>();*/

tests/functional/mhd_harris/harris.py

Lines changed: 49 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,29 +1,57 @@
1+
#!/usr/bin/env python3
2+
import os
3+
14
import numpy as np
2-
import pyphare.mock_mhd_simulator.mhd_model as m
3-
import pyphare.mock_mhd_simulator.simulation as s
4-
from pyphare.mock_mhd_simulator.simulator import MHDMockSimulator
5+
import pyphare.pharein as ph
6+
from pyphare.cpp import cpp_lib
7+
from pyphare.simulator.simulator import Simulator, startMPI
8+
9+
os.environ["PHARE_SCOPE_TIMING"] = "1" # turn on scope timing
10+
11+
ph.NO_GUI()
12+
cpp = cpp_lib()
13+
startMPI()
14+
15+
diag_outputs = "phare_outputs/high"
16+
final_time = 10
17+
time_step = 0.0005
18+
time_step_nbr = int(final_time / time_step)
19+
20+
dumpfrequency = 2000
21+
dt = dumpfrequency * time_step
22+
timestamps = dt * np.arange(int(time_step_nbr / dumpfrequency) + 1)
523

624

725
def config():
8-
cells = (250, 250)
9-
dl = (0.1, 0.1)
10-
11-
sim = s.Simulation(
12-
ndim=2,
13-
order=2,
14-
timestep=0.002,
15-
final_time=10,
26+
cells = (500, 500)
27+
dl = (0.05, 0.05)
28+
29+
sim = ph.Simulation(
30+
smallest_patch_size=15,
31+
# largest_patch_size=25,
32+
time_step_nbr=time_step_nbr,
33+
time_step=time_step,
1634
cells=cells,
1735
dl=dl,
18-
origin=(0.0, 0.0),
36+
refinement="tagging",
37+
max_mhd_level=1,
38+
max_nbr_levels=1,
39+
hyper_resistivity=0.0,
40+
resistivity=0.0,
41+
diag_options={
42+
"format": "phareh5",
43+
"options": {"dir": diag_outputs, "mode": "overwrite"},
44+
},
45+
strict=True,
1946
eta=0.0,
2047
nu=0.0,
2148
gamma=5.0 / 3.0,
22-
reconstruction="wenoz",
49+
reconstruction="weno3",
2350
limiter="",
2451
riemann="rusanov",
25-
time_integrator="tvdrk3",
52+
mhd_timestepper="tvdrk3",
2653
hall=True,
54+
model_options=["MHDModel"],
2755
)
2856

2957
Lx = cells[0] * dl[0]
@@ -83,13 +111,18 @@ def bz(x, y):
83111
def p(x, y):
84112
return 1.0 - (bx(x, y) ** 2 + by(x, y) ** 2) / 2.0
85113

86-
m.MHDModel(density=density, vx=vx, vy=vy, vz=vz, bx=bx, by=by, bz=bz, p=p)
114+
ph.MHDModel(density=density, vx=vx, vy=vy, vz=vz, bx=bx, by=by, bz=bz, p=p)
115+
116+
ph.ElectromagDiagnostics(quantity="B", write_timestamps=timestamps)
117+
118+
for quantity in ["rho", "V", "P"]:
119+
ph.MHDDiagnostics(quantity=quantity, write_timestamps=timestamps)
87120

88121
return sim
89122

90123

91124
def main():
92-
MHDMockSimulator(config()).run("hall_harris.h5", dumpfrequency=500)
125+
Simulator(config()).run()
93126

94127

95128
if __name__ == "__main__":

0 commit comments

Comments
 (0)