Skip to content

Commit 76d653b

Browse files
UCaromelUCaromel
committed
solved some problem for J in 1d with multistep timeintegrators, as well as a timeintegrator refactor error
1 parent 3643435 commit 76d653b

File tree

8 files changed

+65
-51
lines changed

8 files changed

+65
-51
lines changed

‎src/amr/messengers/mhd_messenger.hpp‎

Lines changed: 9 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -353,12 +353,6 @@ namespace amr
353353

354354
bool isRegriddingL0 = levelNumber == 0 and oldLevel;
355355

356-
for (auto& patch : *level)
357-
{
358-
auto _ = resourcesManager_->setOnPatch(*patch, mhdModel.state.J);
359-
mhdModel.state.J.zero();
360-
}
361-
362356
magneticRegriding_(hierarchy, level, oldLevel, initDataTime);
363357
densityInitRefiners_.regrid(hierarchy, levelNumber, oldLevel, initDataTime);
364358
momentumInitRefiners_.regrid(hierarchy, levelNumber, oldLevel, initDataTime);
@@ -394,12 +388,6 @@ namespace amr
394388
densityInitRefiners_.fill(levelNumber, initDataTime);
395389
momentumInitRefiners_.fill(levelNumber, initDataTime);
396390
totalEnergyInitRefiners_.fill(levelNumber, initDataTime);
397-
398-
for (auto& patch : level)
399-
{
400-
auto _ = resourcesManager_->setOnPatch(*patch, mhdModel.state.J);
401-
mhdModel.state.J.zero();
402-
}
403391
}
404392

405393
void firstStep(IPhysicalModel& model, SAMRAI::hier::PatchLevel& level,
@@ -518,13 +506,19 @@ namespace amr
518506
// likely
519507
void registerGhostComms_(std::unique_ptr<MHDMessengerInfo> const& info)
520508
{
509+
// static refinement for J and E because in MHD they are temporaries, so keeping there
510+
// state updated after each regrid is not a priority. However if we do not correctly
511+
// refine on regrid, the post regrid state is not up to date (in our case it will be nan
512+
// since we nan-initialise) and thus is is better to rely on static refinement, which
513+
// uses the state after computation of ampere or CT.
521514
elecGhostsRefiners_.addStaticRefiners(info->ghostElectric, EfieldRefineOp_,
522515
info->ghostElectric,
523516
nonOverwriteInteriorTFfillPattern);
524517

525-
currentGhostsRefiners_.addTimeRefiners(info->ghostCurrent, info->modelCurrent,
526-
Jold_.name(), EfieldRefineOp_, vecFieldTimeOp_,
527-
nonOverwriteInteriorTFfillPattern);
518+
currentGhostsRefiners_.addStaticRefiners(info->ghostCurrent, EfieldRefineOp_,
519+
info->ghostCurrent,
520+
nonOverwriteInteriorTFfillPattern);
521+
528522

529523
rhoGhostsRefiners_.addTimeRefiners(info->ghostDensity, info->modelDensity,
530524
rhoOld_.name(), mhdFieldRefineOp_, fieldTimeOp_,

‎src/amr/messengers/mhd_messenger_info.hpp‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
#ifndef PHARE_MHD_MESSENGER_INFO_HPP
22
#define PHARE_MHD_MESSENGER_INFO_HPP
33

4-
#include "core/data/vecfield/vecfield.hpp"
54
#include "core/numerics/godunov_fluxes/godunov_utils.hpp"
65
#include "messenger_info.hpp"
76

‎src/amr/physical_models/mhd_model.hpp‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -107,10 +107,6 @@ void MHDModel<GridLayoutT, VecFieldT, AMR_Types, Grid_t>::initialize(level_t& le
107107
auto _ = this->resourcesManager->setOnPatch(*patch, state);
108108

109109
state.initialize(layout);
110-
// data initialized to NaN on construction
111-
// and in 1D Jx is not worked on in Ampere so
112-
// we need to zero J before anything happens
113-
state.J.zero();
114110
}
115111
resourcesManager->registerForRestarts(*this);
116112
}

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

Lines changed: 8 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -87,35 +87,33 @@ class BaseMHDTimestepper
8787
}
8888
}
8989

90-
void accumulateButcherFluxes_(MHDModel& model, auto& fluxes, auto& level)
90+
void accumulateButcherFluxes_(MHDModel& model, auto& E, auto& fluxes, auto& level,
91+
double const coef = 1.0)
9192
{
9293
for (auto& patch : level)
9394
{
9495
auto const& layout = amr::layoutFromPatch<GridLayoutT>(*patch);
95-
auto _ = model.resourcesManager->setOnPatch(*patch, butcherFluxes_, butcherE_, fluxes,
96-
model.state.E);
96+
auto _
97+
= model.resourcesManager->setOnPatch(*patch, butcherFluxes_, butcherE_, fluxes, E);
9798

9899
evalFluxesOnGhostBox(
99100
layout,
100101
[&](auto& left, auto const& right, auto const&... args) mutable {
101-
left(args...) += right(args...);
102+
left(args...) += right(args...) * coef;
102103
},
103104
butcherFluxes_, fluxes);
104105

105106

106107
layout.evalOnGhostBox(butcherE_(core::Component::X), [&](auto const&... args) mutable {
107-
butcherE_(core::Component::X)(args...)
108-
+= model.state.E(core::Component::X)(args...);
108+
butcherE_(core::Component::X)(args...) += E(core::Component::X)(args...) * coef;
109109
});
110110

111111
layout.evalOnGhostBox(butcherE_(core::Component::Y), [&](auto const&... args) mutable {
112-
butcherE_(core::Component::Y)(args...)
113-
+= model.state.E(core::Component::Y)(args...);
112+
butcherE_(core::Component::Y)(args...) += E(core::Component::Y)(args...) * coef;
114113
});
115114

116115
layout.evalOnGhostBox(butcherE_(core::Component::Z), [&](auto const&... args) mutable {
117-
butcherE_(core::Component::Z)(args...)
118-
+= model.state.E(core::Component::Z)(args...);
116+
butcherE_(core::Component::Z)(args...) += E(core::Component::Z)(args...) * coef;
119117
});
120118
}
121119
}

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ class EulerIntegrator : public BaseMHDTimestepper<MHDModel>
2929

3030
euler_(model, state, state, fluxes, bc, level, currentTime, newTime);
3131

32-
this->accumulateButcherFluxes_(model, fluxes, level);
32+
this->accumulateButcherFluxes_(model, state.E, fluxes, level);
3333
}
3434

3535
using Super::allocate;

‎src/core/numerics/MHD_equations/MHD_equations.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -160,7 +160,7 @@ class MHDEquations
160160
F_B.z += -J.x * coef;
161161
F_Etot += (J.z * B.x - J.x * B.z) * coef;
162162
}
163-
if constexpr (direction == Direction::Y)
163+
if constexpr (direction == Direction::Z)
164164
{
165165
F_B.x += -J.y * coef;
166166
F_B.y += J.x * coef;

‎src/core/numerics/ampere/ampere.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,9 @@ class Ampere_ref
6767
{
6868
auto const& [_, By, Bz] = B();
6969

70+
if constexpr (dimension == 1)
71+
Jx(ijk...) = 0.0;
72+
7073
if constexpr (dimension == 2)
7174
Jx(ijk...) = layout_.template deriv<Direction::Y>(Bz, {ijk...});
7275

‎src/python3/cpp_mhd_parameters.hpp‎

Lines changed: 43 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -208,28 +208,28 @@ constexpr void declare_all_mhd_params(py::module& m)
208208

209209
std::string variant_name = "euler_constant_rusanov";
210210
std::string full_type = type_name + "_" + variant_name;
211-
212-
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
213-
ReconstructionType::Constant, SlopeLimiterType::count,
214-
RiemannSolverType::Rusanov, false, false, false>::declare_sim(m, full_type);
215-
216-
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
217-
ReconstructionType::Constant, SlopeLimiterType::count,
218-
RiemannSolverType::Rusanov, false, false, false>::declare_etc(m, full_type);
219-
220-
// variant_name = "euler_constant_rusanov_hall";
221-
// full_type = type_name + "_" + variant_name;
222211
//
223212
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
224213
// ReconstructionType::Constant, SlopeLimiterType::count,
225-
// RiemannSolverType::Rusanov, true, false, false>::declare_sim(m,
214+
// RiemannSolverType::Rusanov, false, false, false>::declare_sim(m,
226215
// full_type);
227216
//
228217
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
229218
// ReconstructionType::Constant, SlopeLimiterType::count,
230-
// RiemannSolverType::Rusanov, true, false, false>::declare_etc(m,
219+
// RiemannSolverType::Rusanov, false, false, false>::declare_etc(m,
231220
// full_type);
232221

222+
variant_name = "euler_constant_rusanov_hall";
223+
full_type = type_name + "_" + variant_name;
224+
225+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
226+
ReconstructionType::Constant, SlopeLimiterType::count,
227+
RiemannSolverType::Rusanov, true, false, false>::declare_sim(m, full_type);
228+
229+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::Euler,
230+
ReconstructionType::Constant, SlopeLimiterType::count,
231+
RiemannSolverType::Rusanov, true, false, false>::declare_etc(m, full_type);
232+
233233
// variant_name = "tvdrk3_wenoz_rusanov";
234234
// full_type = type_name + "_" + variant_name;
235235
//
@@ -242,7 +242,7 @@ constexpr void declare_all_mhd_params(py::module& m)
242242
// ReconstructionType::WENOZ, SlopeLimiterType::count,
243243
// RiemannSolverType::Rusanov, false, false, false>::declare_etc(m,
244244
// full_type);
245-
//
245+
246246
// variant_name = "tvdrk3_wenoz_rusanov_hall";
247247
// full_type = type_name + "_" + variant_name;
248248
//
@@ -268,17 +268,41 @@ constexpr void declare_all_mhd_params(py::module& m)
268268
// ReconstructionType::Linear, SlopeLimiterType::VanLeer,
269269
// RiemannSolverType::Rusanov, false, false, false>::declare_etc(m,
270270
// full_type);
271+
272+
variant_name = "tvdrk2_linear_vanleer_rusanov_hall";
273+
full_type = type_name + "_" + variant_name;
274+
275+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK2,
276+
ReconstructionType::Linear, SlopeLimiterType::VanLeer,
277+
RiemannSolverType::Rusanov, true, false, false>::declare_sim(m, full_type);
278+
279+
RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK2,
280+
ReconstructionType::Linear, SlopeLimiterType::VanLeer,
281+
RiemannSolverType::Rusanov, true, false, false>::declare_etc(m, full_type);
282+
283+
// variant_name = "tvdrk3_weno3_rusanov";
284+
// full_type = type_name + "_" + variant_name;
285+
//
286+
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
287+
// ReconstructionType::WENO3, SlopeLimiterType::count,
288+
// RiemannSolverType::Rusanov, false, false, false>::declare_sim(m,
289+
// full_type);
271290
//
272-
// variant_name = "tvdrk2_linear_vanleer_rusanov_hall";
291+
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
292+
// ReconstructionType::WENO3, SlopeLimiterType::count,
293+
// RiemannSolverType::Rusanov, false, false, false>::declare_etc(m,
294+
// full_type);
295+
296+
// variant_name = "tvdrk3_weno3_rusanov_hall";
273297
// full_type = type_name + "_" + variant_name;
274298
//
275-
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK2,
276-
// ReconstructionType::Linear, SlopeLimiterType::VanLeer,
299+
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
300+
// ReconstructionType::WENO3, SlopeLimiterType::count,
277301
// RiemannSolverType::Rusanov, true, false, false>::declare_sim(m,
278302
// full_type);
279303
//
280-
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK2,
281-
// ReconstructionType::Linear, SlopeLimiterType::VanLeer,
304+
// RegistererSelector<Dimension, InterpOrder, NbRefinedParts, TimeIntegratorType::TVDRK3,
305+
// ReconstructionType::WENO3, SlopeLimiterType::count,
282306
// RiemannSolverType::Rusanov, true, false, false>::declare_etc(m,
283307
// full_type);
284308

0 commit comments

Comments
 (0)